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ABSTRACT 


Analysis  of  the  transient  time  response  of  many  modern  engineering 
systems  requires  simulation  by  means  of  direct  time  integration  of  sets  of 
coupled-field  equations.  In  the  past  the  simulation  of  such  coupled 
problems  has  created  difficulties  which  are  not  seen  in  analogous  single¬ 
field  problems.  Perhaps  the  most  important  of  these  difficulties  is  the  fact 
that  existing  engineering  analysis  packages  have  been  oriented  to  solving 
only  single-field  problems.  This  suggests  one  particularly  attractive 
option:  the  use  of  existing  single -field  analysis  software  while  imposing 
certain  modular  requirements.  If  modularity  can  be  maintained,  direct 
integration  of  the  equations  can  proceed  by  sequential  or  parallel 
operation  of  the  separate  analyzers. 

This  research  develops  the  theory  by  which  the  numerical 
simulations  of  two  or  more  separate  fields  may  be  combined  to  solve  the 
coupled-field  problem.  This  theory  allows  simulations  to  be  used  with 
little  or  no  change  by  considering  the  constraints  that  provide  for  field 
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coupling.  The  development  takes  into  account  the  various  numerical 
methods  which  may  be  used  by  individual  simulations  to  solve  their 
separate  problems.  Above  all  this  paper  does  not  seek  to  suggest  that 
simulation  coupling  is  the  best  or  foremost  simulation  methodology 
available,  merely  that  it  is  a  viable  and  cost  saving  alternative  to  solving 
the  larger,  more  involved  coupled-field  problem. 

The  specific  discipline  of  focus  will  be  multibody  dynamics.  The  goal 
will  be  to  show  that  it  is  possible  to  couple  multiple  single  body  dynamic 
analysis  packages  and  come  up  with  solutions  comparable  to  the  full 
multibody  dynamic  case.  The  advantages  of  such  a  methodology  are 
readily  seen  here.  The  single  body  equations  of  motion  are  fairly  simple 
and  a  great  deal  of  software  exists  to  produce  transient  analyses  of  such 
bodies.  This  is  compared  to  the  multibody  analysis  where  the  equations  of 
motion  are  not  readily  derived,  few  reliable  analyzers  exist,  and  changes  in 
the  model  may  require  a  great  deal  of  change  to  the  analysis  package. 

The  paper  is  organized  into  three  main  sections.  The  first  section 
deals  with  the  tools  necessary  to  evaluate  the  simulation  coupling  method. 
The  second  part  introduces  this  method  as  well  as  dealing  with  stability 
and  accuracy  of  the  algorithm  presented.  The  final  section  of  the  paper 
involves  numerical  examples  which  demonstrate  the  strengths  and 
weaknesses  of  this  theory. 

Thesis  Supervisor:  Dr.  David  S.  Kang 

Technical  Staff,  Charles  Stark  Draper  Laboratory 
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Chapter  One 


Introduction 


1.1  Background 

Flexible  multibody  problems  have  become  increasingly 
important  in  recent  years,  primarily  to  support  the  design  and 
deployment  of  large  space  structures.  The  numerical  simulation  of 
these  problems  becomes  especially  important  in  space  applications 
where  extensive  laboratory  research  is  impractical.  Flexible 
multibody  dynamics,  however,  is  only  a  subset  of  the  larger  group  of 
coupled-field  problems  whose  numerical  solution  has  been  the 
subject  of  a  great  deal  of  study  over  the  past  decade.  Coupled-field 
systems  are  readily  solved  single-field  subsystems  linked  together 
by  constraints  or  interface  boundaries;  systems  whose  solution  is 
made  difficult  by  the  size  of  the  combined  problems,  the  widely 
varying  time  response  characteristics  that  the  combined  subsystems 
may  have,  and  the  fact  that  most  analysis  software  currently 
available  is  written  for  single-field  analysis. 

Traditional  solutions  to  the  coupled-field  problem  involve 
solving  the  system  as  a  whole,  either  through  field  elimination  or  by 
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simultaneous  solution.  Recently,  partitioned  solution  methods  have 
received  a  great  deal  consideration  for  their  ability  to  divide  the 

problem  into  pieces  which  may  be  dealt  with  individually.  A  more 
extreme  form  of  partitioning  exists  in  the  form  of  simulation 

coupling,  where  the  constraint  or  boundary  interface  effects  are 

combined  with  single  field  analyzers  to  solve  coupled  systems. 

The  major  advantage  of  the  simulation  coupling  method  is  its 
maximal  use  of  available  software  and  existing  analyses.  This  is 
especially  true  for  coupled  structural  problems.  Transient  analysis  of 
such  problems  are  of  particular  interest  in  the  design  process  when 
variations  of  secondary  structural  systems  need  to  be  coupled  to  a 
primary  structure.  Analysis  of  the  loading  due  to  different  satellites 
carried  in  the  space  shuttle's  main  bay  and  the  construction  and 

deployment  of  laboratory  and  habitation  modules  connected  to  the 
space  station's  main  structure  are  two  examples  where,  if  certain 
stability  and  accuracy  needs  can  be  reached,  there  would  be  an 
advantage  to  coupling  existing  simulations  instead  of  using  the  more 
traditional  methods  to  re-solve  the  problem  for  each  variation. 

1.2  Motivation  for  Current  Work 

Analysis  of  the  transient  time  response  of  many  modern 
engineering  systems  requires  simulation  by  means  of  direct  time 
integration  of  sets  of  coupled-field  equations.  The  simulation  of  such 
coupled  problems  has  created  difficulties  which  are  not  normally 
seen  in  analogous  single-field  problems.  Perhaps  the  most  important 
of  these  difficulties  is  the  fact  that  existing  engineering  analysis 
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packages  have  been  oriented  to  solving  single-field  problems.  This 
suggests  one  particularly  attractive  answer:  the  use  existing  single¬ 
field  analysis  software  while  imposing  certain  modular  requirements. 
If  modularity  can  be  maintained,  direct  integration  of  the  equations 
can  proceed  by  sequential  or  parallel  operation  of  separate  analyzers. 

The  motivation  for  this  thesis  is  the  development  of  a  theory 
by  which  the  numerical  simulations  of  two  or  more  separate  fields 
may  be  combined  to  solve  the  coupled-field  problem.  This  theory 
should  allow  simulations  to  be  used  without  change  through 
consideration  of  the  constraints  that  provide  for  the  field  coupling. 
This  theory  should  account  for  the  various  numerical  methods  which 
may  be  used  by  the  individual  simulations  to  solve  their  separate 
problems.  It  should  be  noted  that  this  thesis  is  aimed  at  solving 
coupled  problems  where  at  least  one  field  is  structural  in  nature. 
Even  then  the  arguments  and  theory  contained  within  this  thesis  are 
aimed  at  the  more  physical  interaction  problems  (ie.  fluid-structure 
and  structure-structure  problems  versus  control-structure  and 
thermal-mechanical  problems). 

1.3  Restrictions  Considered 

In  developing  the  theory  described  earlier  certain  restrictions 
of  the  general  problem  are  assumed.  First,  the  constraints  which 
couple  the  separate  fields  are  assumed  to  exist  as  two-way 
interactions  between  two  individual  fields.  For  this  reason  the 
coupled  problem  can  be  examined  by  considering  only  two  fields. 
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Additional  fields  can  be  dealt  with  in  a  similar  fashion.  Figure  1.3-1 
shows  the  allowable  interactions  for  a  system  of  three  coupled  fields. 

Separate  Fields 


Figure  1.3-1  Coupled  System  Interaction 

It  is  also  assumed  that  the  equations  of  motion  for  the  separate 
fields  are  initially  describable  in  a  semi-discrete  second  order  form: 

Mx  +  Cx  +  Kx  =  R(t) 

This  is  not  an  overly  limiting  assumption  since  most  fields  in  coupled 
problems  are  easily  placed  in  this  form. 

1.4  Overview 

The  remainder  of  this  thesis  is  organized  in  the  following 
manner.  Chapter  two  reviews  the  available  methods  of  solving 
coupled  field  problems  along  with  the  benefits  and  shortcomings  of 
each  method.  Field-elimination  and  simultaneous  solutions  are 
discussed  as  well  as  partitioned  solution  methods,  where  simulation 
coupling  is  introduced  as  a  special  case  of  these  methods. 
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Chapter  Three  describes  some  of  the  tools  essential  to  the 
development  of  the  theory  for  coupling  simulations,  namely 
integration  procedures  and  stability  analysis  procedures. 

Chapter  Four  sets  forth  the  simulation  coupling  procedures  in 
addition  to  discussing  how  to  deal  with  the  constraints  which  couple 
the  separate  fields.  The  stability  and  accuracy  of  the  coupling 
procedures  are  also  discussed. 

Chapters  Five  and  Six  set  forth  numerical  examples  designed  to 
illuminate  points  made  about  the  simulation  coupling  theory 
developed  in  chapter  five.  Chapter  Five  involves  small  scale 
examples  designed  to  highlight  specifics  on  the  stability  and  accuracy 
of  this  method.  Chapter  Six  is  dedicated  to  the  large  scale  example, 
namely,  coupling  single  body  simulations  to  form  the  assembly 
complete  form  of  the  space  station;  a  problem  at  the  scale  to  which 
simulation  coupling  is  specifically  directed. 

Chapter  Seven  concludes  this  thesis  by  reviewing  work  done, 
major  findings  included  within,  as  well  as  outlining  future  work  for 
consideration  in  this  area. 
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Chapter  Two 


Solution  Methods  of 
Coupled-Field  Problems 


This  chapter  presents  a  basic  overview  of  the  primary  means 
by  which  coupled-field  problems  are  solved.  These  problems  are 
often  sufficiently  large  and  complex  that  the  only  feasible  solution 
process  is  direct  integration,  which  leaves  the  major  question  of  what 
form  to  place  the  equations  of  motion  in  order  to  carry  out  the 

integration.  Field  elimination  and  simultaneous  solution  are  standard 
methods  by  which  these  problems  are  currently  handled  although 

both  have  serious  drawbacks  which  hinder  their  performance.  More 
recently  partitioned  solutions  have  been  suggested  as  a  possible 
alternative  [2.1-2.10].  A  large  number  of  partitioning  schemes  have 
been  recommended  in  literature  and  some  of  the  more  popular 

methods  will  be  introduced  here.  At  the  end  of  this  chapter 
simulation  coupling  is  shown  as  a  natural  extension  of  partitioned 

solutions. 
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2.1  Idealization  of  the  Coupled-Field  Problem 

Consider  an  arbitrary  domain,  S,  shown  in  Figure  2.1-1. 
Allowing  that  S  is  the  domain  of  a  two-field  coupled  problem  it  may 
be  decomposed  into  three  distinct  subdomains:  Sx  and  Sy,  the 
separate  single  field  subdomains;  and  Si,  the  interface  subdomain. 


If  the  effects  of  Si  are  ignored  then  the  single  fields  may  be 
modelled  by  finite  difference  or  finite  element  methods  as 
semidiscrete,  linear,  second  order  matrix  differential  equations 

Mxx  +  Cxx  +  Kxx  =  Rx 

Myy +  Cyy  +  Kyy  =  Ry  (2.1-1) 

where 

x  =  x(t),y  =  y(t)  :  Separate  field  response  vectors 
M  x,  M  y  :  General  mass  matrix 
Cx,  Cy  :  General  damping  matrix 
Kx,  Ky  :  General  stiffness  matrix 
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Rx  =  Rx(t),Ry  =  R y(t)  :  Separate  field  vector  forcing  functions 
and  ( * )  denotes  time  integration 

The  fields  expressed  in  equation  (2.1-1)  are  coupled  by  the  physical 
constraints  and  interface  conditions  of  Si  which  are  expressed  as 

®(x,y,x,y)  =  0  (2.1-2) 

Note  that  the  constraints  which  lead  to  coupled-field  problems  are 
generally  not  assumed  to  be  time  varying  functions. 

Through  appropriate  use  of  Lagrange  multipliers,  penalty 
formulations,  or  other  methods  equations  (2.1-1)  and  (2.1-2)  may  be 
combined  as 


Mxx  +  Cxx  +  Kxx  =  Rx  +  fcx 

Myy  +  Cyy  +  Kyy  =  Ry  +  fcy  (2.1-3) 

where  fcx  and  fcy  are  effective  constraint  forces  used  to  correct  the 
separate  field  response. 

In  the  case  where  the  interactions  are  linear  in  nature  then 
(2.1-3)  can  be  rearranged  as  follows: 

Mxxx  +  Cxxx  +  Kxxx  =  Rx  -  CXyy  -  KXyy 

Myyy  +  Cyyy  +  Kyyy  =  Ry  ■  CyXx  -  KyXx  (2.1  -4) 

or  if  u  =  [  xT  yT  ]T 

Mil  +  Cii  +  Ku  =  R  (2.1-5) 


where 
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Mxx  0  Cxx  CXy  Kxx  Kxy 

M=  C  =  K - 

[  0  MyyJ’  ^  LCyxCyyJ’  ^  [  Kyx  Kyy  J  ’ 

andR=[RxT  RyTlT 

Note:  Kx  may  be  modified  by  terms  from  fcx  in  (2.1-3)  in  order  to 
obtain  Kxx  in  (2.1-4).  Similar  effects  may  be  seen  in  Kyy,  Cxx,  Cyy. 
However,  Mxx  =  Mx  and  Myy  =  My. 

2.2  Field  Elimination 

Field  elimination  is  aimed  at  the  elimination  of  the  interaction 
terms  through  substitutions  from  the  other  available  equations.  This 
method  is  usually  suggested  when  the  response  of  one  field  is 
considered  to  be  more  important  than  the  other.  This  process 
reduces  the  number  of  states  associated  with  each  problem  to  only 
one  field  but  has  several  considerable  drawbacks  which  should  be 
highlighted  in  the  following  example. 

Consider  the  simpler  system  of  the  form  of  (2.1-4) 

Mxxx  +  Cxx*  +  KXxX  =  Rx  -  Kxyy 

Myyjr  +  Cyyy  +Kyyy  =  Ry  -  Kyxx  (2.2-1) 

Transforming  these  equations  by  means  of  the  Laplace  variable 

(M„s2  +  C„s  +  K JX(s)  =  Ms)  -  KxyY(s) 

(MyyS2  +  CyyS  +  Kyy)Y(s)  =  Ry(s)  -  KyXX(s)  (2.2-2) 

Eliminating  Y(s)  from  (2.2-2a)  using  (2.2-2b) 
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[(lUyyS2  +  CyyS  +  Ky  yjx  xy ' '(iVI  XXS2  +  C„S  +  Kxx)  -  Kyjx(s)  = 

(MyyS2  +  CyyS  +  K  y  yjk  *  y’ ‘  R  ,  (s)  -  Ry(s)  (2-2'3) 

Multiplying  through 

{(MyyKXy-1MxJs4+(CyyKXy-1MXX  +  MyyKXy-1CXX)s3 

+  (MyyKXy  KXX  +  CyyKXy  CXX  +  KyyKXy  Mxx)s2 

+  (KyyKxy-ICXX+CyyKXy-IKxx)s+(KyyKXy-1Kxx-KyX)}X(s)  (2.2-4) 

=  (MyyS2  +  CyyS  +  Ryyjx xy  RX(S)  -  Ry(s) 

Use  of  the  inverse  Laplace  transform  returns  a  differential 
expression 

(MyyKXy  Mxx)  X  +(CyyKXy  MXX  +  MyyKXy  CXX)j[ 

+  (MyyKXy  KXX+CyyKXy  CXX+KyyKXy  MXX)jC 
+  (KyyKXy-‘Cxx  +  C  y  y  Kxy“ '  K  X  X)  X  +  (x  yy  K  X  y“ 1  K  X  x  -  KyX)x  (2'2'5) 

=  ^®yy^xy  !Rx  +  CyyKXy  RX  +  KyyKXy  Rx  -  Ry 

An  expression  similar  to  (2.2-4)  can  be  found  for  the  variable  y(t). 
The  following  disadvantages  can  be  seen  in  this  example: 


Higher-order  derivatives  are  introduced  for  which  there  are  no 
readily  available  integrators. 

-  Additional  initial  conditions  are  required. 

Additional  derivatives  of  the  forcing  functions  are  required. 

The  sparsity  (and  possibly  symmetry)  of  the  attendant  matrices  is 
lost. 

There  is  little  possible  use  of  existing  software  or  single  field 
analyzers 
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The  above  disadvantages  become  more  serious  as  the  complexity 
of  the  problem  increases. 

Although  field  elimination  does  succeed  in  reducing  the  states 
associated  with  the  problem  it  does  so  at  the  cost  of  introducing 
further  difficulties  not  easily  solved  at  this  time.  Furthermore,  every 
time  a  secondary  field  to  be  coupled  with  the  primary  changes  the 
process  must  be  repeated  to  solve  the  new  problem. 

2.3  Simultaneous  Solutions 

In  simultaneous  solution  methods  the  coupled  equations  of 
motion  are  integrated  as  a  single  large  second-order  problem  in  the 
form  of  equation  (2.1-5).  By  integrating  all  the  equations  in  second- 
order  form,  as  opposed  to  field  elimination  procedures,  higher  order 
derivatives  and  additional  initial  conditions  are  not  required.  In 
addition  there  are  accepted  integration  methods  for  use  on  second- 
order  differential  systems.  However,  to  be  placed  in  the  form  of 
(2.1-5)  the  interactions  must  be  linear  and  there  is  still  no  possibility 
of  using  existing  single  field  analyzers. 

Another  shortcoming  of  this  method  is  the  requirement  of 
treatment  in  fully  explicit  or  fully  implicit  form.  This  specification 
carries  with  it  the  following  problems.  Since  coupled  problems 
typically  have  very  diverse  time  characteristics,  the  stability  limits 
on  the  step  size  tend  to  be  unreasonably  restrictive.  This  is 
especially  true  if  the  problem  includes  rigid  effects  or  incorporates 
the  constraints  by  penalty  formulations.  Also  the  interaction  tends 


30 


CHAPTER  TWO:  SOLUTION  METHODS  OF  COUPLED-FIELD  PROBLEMS 


to  produce  extremely  large  bandwidths  in  the  associated  coefficient 
matrices.  As  a  result  the  solution  of  realistic  three  dimensional 
problems  becomes  rapidly  prohibitive  due  to  the  number  of 
calculations  necessary  to  set  up  and  solve  equations  involving  these 
matrices. 

So  despite  the  superiority  of  simultaneous  solutions  to  field 
elimination  there  still  exists  the  problem  of  setting  up  and  solving 
large  set  of  coupled  equations,  along  with  a  continued  lack  of 
modularity. 

2.4  Partitioned  Solutions 

Partitioned  solution  methods  are  based  on  dividing  the  system 
matrices  of  equation  (2.1-5)  into  two  parts, 

K  =  Ki  +  KE  and  C  =  O  +  CE  (2.4-1) 

where  K1  and  C1  are  the  implicit  portions  of  the  partition  and  KE  and 
CE  are  the  explicit  portions.  It  should  be  mentioned  that  the  entire 
mass  matrix  must  be  contained  in  the  implicit  portion  of  the  partition 
(see  [2.6]).  In  a  partitioned  solution  the  explicit  portion,  combined 
with  a  predictor,  acts  like  an  applied  force  input  to  the  differential 
equation.  There  are  two  things  which  define  a  partition  method:  the 
partitioning  strategy  and  the  particular  partition  used  to  divide  the 
system  matrices. 

The  partitioning  strategy  is  defined  by  the  point  at  which  the 
partitioning  occurs.  At  some  point  in  order  to  solve  the  differential 
system  a  numerical  integration  scheme  must  be  applied.  If  the 
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scheme  is  applied  and  then  the  resulting  matrix  equation  is 
partitioned  the  strategy  is  call  algebraic  partitioning.  If  the  partition 
is  applied  before  the  integration  scheme,  then  a  differential 
partitioning  strategy  results.  The  work  of  Felippa  and  Park  [2.6] 
shows  a  strong  connection  between  the  stability  of  the  method  and 
the  partitioning  strategy  used. 

For  a  two-field  system  there  are  sixteen  possible  ways  to  fully 
partition  the  system  matrices  and  six  of  these  simply  represent  field 
switches.  The  ten  remaining  unique  partitions  for  the  explicit 
portion,  KE  are 

‘ro  o]  2r  o  o  i  3r  o  o  1 4r  o  o  1 5r  o  Kxy ' 

-0  oJ  L  0  KyyJ  [ Kyx  0  L  Kyx  Kyy  J  [  0  Kyy  J(2.4-2) 

6  Kxx  0  7  0  Kxy  *  Kxx  Kxy  ^  Kxx  0  Kxx  Kxy 

0  Kyy  Kyx  0  KyX  0  KyX  Kyy  KyX  Kyy 

The  first  and  last  cases  correspond  to  the  limiting  fully  implicit  and 
fully  explicit  simultaneous  solutions,  respectively.  Partitions 
numbers  2,  4,  5,  8,  and  9  are  all  implicit-explicit  partitions  with  2 
and  4  being  more  widely  used,  as  will  be  discussed  later.  Also 

number  3  is  of  particular  importance;  it  is  referred  to  as  a  staggered 
partition  and  warrants  further  consideration.  A  more  complete 
analysis  of  all  the  available  partitions  is  contained  in  [2.7]. 

Implicit-explicit  partitions  were  first  suggested  for  use  in 

structural  dynamics  by  Belytschko  and  Mullen  [2. 1,2.2]  in  problems 
where  the  mesh  had  two  distinct  set  of  time  characteristics.  The 
specific  application  in  mind  were  fluid-structure  problems  where  a 
very  large,  slow  responding  mesh  (fluid)  was  coupled  with  a  smaller, 
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quicker  responding  mesh  (structure).  Since  the  fluid  mesh  is  very 
large  an  explicit  method  was  desired  for  improved  computational 
efficiency,  but  the  large  range  of  response  frequencies  in  the 
structure  called  for  an  implicit  method  for  stability.  By  introducing  a 
boundary  field,  thus  making  it  a  three-field  problem,  and  a  three- 
field  partition  equivalent  to  4  the  fluid  was  dealt  with  explicitly  and 
the  structure  implicitly.  This  form  of  partition  is  known  as  node-by- 
node  implicit-explicit  (IE)  partitioning. 

Shortly  thereafter  Hughes  and  Liu  [2. 3,2. 4]  introduced  the 
three-field  partition  based  on  2.  This  model  defined  the  elements  as 
either  implicit  or  explicit.  The  element-by-element  IE  partitioning  is 
easier  to  implement  but  may  be  more  computationally  intensive  for 
large  boundary  problems.  However  both  element-by-element  and 
node-by-node  IE  partitioning  still  retain  little  modularity  in  solving 
the  problem. 

Perhaps  the  most  extensively  used  partition  is  the  staggered 
partition.  It  has  thus  far  been  effectively  applied  to  fluid-structure 
interaction  problems  [2.10],  control-structure  interaction  problems 
[2.9]  and  multibody  dynamic  simulation  [2.8].  Staggered  solutions 
predict  one  field  (in  this  case  x)  and  use  the  prediction  to  solve  the 
second  field  (y).  The  first  field  is  then  solved  using  the  solution  of  the 
second  in  an  implicit  fashion.  This  process  is  depicted  below  in 
Figure  2.4-1  where  EP  is  an  explicit  prediction  flow  and  IS  is  an 
implicit  solution  flow.  Although  it  has  been  successful  in  allowing 
more  use  of  single  field  analysis  software  [see  2.6],  it  still  does  not 
come  through  with  the  desired  modularity. 
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2.5  Simulation  Coupling 

Simulation  coupling  does  not  fall  under  one  of  the  ten 
partitions  already  mentioned  since  it  is  not  a  complete  partition. 
Instead  of  starting  from  (2.1-5)  the  partition  is  formed  from  the 
terms  in  (2.1-3) 


Mxx  +  Cxx  +  Kxx  —  Rx  +  fcx 
Myy  +  Cyy  +  Kyy  =  Ry  +  fcy 


If  these  equations  are  manipulated  into  the  form  of  (2.1-5)  they  look 
like 


Mx 

0 


cx 

0 


ii  + 


Kx 

0 


0 

Ky 


u  =  R  + 


fcx 
fcy  . 


(2.5-1) 


where  all  terms  are  as  previously  defined. 

Again  by  approximating  the  interaction  terms  as  linear,  the 
effective  constraint  terms  may  be  written 
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Using  these  terms  the  partition  for  simulation  coupling  is 
defined  as  follows: 


Kx 

0 


0 

Ky 


and  Ke  = 


wfC  i/C 
**xx  "xy 

i/C  I^C 

*^yx  **yy 


(2.5-3) 


and  the  damping  terms  are  similarly  partitioned.  The  advantages  to 
using  such  a  partition  are  that  the  implicit  portion  of  the  partition  is 
exactly  the  single  field  problem  with  interaction  effects  neglected. 
This  is  readily  seen  by  examining  the  left-hand  side  of  (2.5-1),  where 
the  system  matrices  are  completely  uncoupled.  The  interaction 
effects  are  handled  exclusively  in  the  explicit  portion  of  the  partition. 
This  solution  form  allows  existing  single  field  analysis  packages  to  be 
applied  directly  with  the  inclusion  of  constraint  terms  as  applied 
forces  being  the  only  modification  necessary.  As  this  is  a  partitioned 
solution  form,  the  partitioning  strategy,  the  chosen  predictors,  as  well 
as  the  details  of  implementation  have  great  impact  on  overall 
stability  and  performance  of  simulation  coupling.  These  effects  are 
addressed  in  the  remaining  parts  of  this  thesis. 
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Chapter  Three 


Methods  for  the  Analysis  of 
Numerical  Integration 


Since  the  primary  means  of  solving  large  sets  of  differential 
equations  is  direct  time  integration,  it  is  important  to  review  the 
basic  theory  and  notation  involved.  This  chapter  outlines  reduced 
order  forms  as  well  as  operational  notation  used  to  determine  the 
stability  of  a  given  procedure.  Additionally,  details  of  computer 
implementation  such  as  computational  paths  and  choice  of  auxiliary 
vectors  are  covered.  Much  of  this  information  was  first  introduced 
by  Jensen  [3.1]  and  is  covered  in  detail  in  a  series  of  papers  by 
Felippa  and  Park  [3.2,  3.3  and  2.6]. 

3.1  Time  Discretization 

The  equation  for  a  general  space-discretized  structural  system 
described  earlier  is 

Mu  +  Cii  +  Ku  =  R  (3.1-1) 
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This  system  may  be  placed  in  reduced  first  order  form  by 
introducing  an  auxiliary  vector,  v  (see  [3.1]) 

v  =  AMii  +  Bu  (3.1-2) 

where  A  and  B  are  suitably  chosen  n  x  n  matrices.  By  manipulating 
(3.1-1)  ii  may  be  eliminated  in  favor  of  v 

AMii  +  ACu  +  AKu  =  AR  ->  v  +(AC  -  B)ti  +  AKu  =  AR  (3.1-3) 

With  the  equations  of  motion  cast  in  first  order  form  [(3.1-2)  and 
(3.1-3)],  numerical  integration  is  carried  out  by  introducing  a  first 
order,  linear  multistep  integration  approximation  for  the  variables  ii 
and  v.  Given  a  constant  stepsize,  h,  the  form  of  such  approximations 
is 

m  m 

2  «iwn-i  =  PiWn-i 

i=0  i=0  (3.1-4) 

The  aj‘s  and  Pi‘s  are  specific  to  each  approximation  and  frequently 

normalized  so  that  ao  =  1.  Also,  Wk  is  a  generic  vector  with  k 
denoting  the  vector  w(t)  at  time  t=tk. 

Note  1:  A  large  number  of  integration  methods  have  been  left 
out  due  to  the  restriction  to  linear  multistep  methods.  Perhaps  the 
most  frequently  used  of  these  are  the  Runge-Kutta  class.  Although 
popular,  this  class  is  not  feasible  for  large  scale  problems  like 
structural  simulation  due  to  their  multiple  derivative  evaluations  per 
time  step  and  the  difficulty  associated  with  analyzing  the  stability  of 
multiple  evaluation  methods. 
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Note  2:  It  is  assumed  that  the  same  integration  approximation 
is  used  for  both  u  and  v.  This  is  consistent  with  the  basic  principles 
of  simulation  coupling;  in  making  maximum  use  of  existing  software 
a  single  favorite  or  best  available  integration  package  would  be 
applied  to  all  equations  within  the  separate  simulations.  If  this  is  not 
the  case  then  (3.1-4)  may  be  replaced  by 

m  m 

X  aiun-i  =  hX 

i=0  i=0 

m  m 

X  aIVn-i  =  hX  Pi^n-i  (3.1-5) 

i=0  i=0 

and  the  separate  terms  may  be  carried  throughout  the  remaining 
equations. 

Removing  the  current  state  term  (t  =  tn)  from  the  past  terms, 
(3.1-4)  may  be  recast 

wn  =  hp0Wn  +  hn  (3.1-6) 

where 

W  W  w  r  .  .  -i 

hn  =  hbn  -  Un  =  h[w  n-l  •  •  •  W  n-mj1 

hPo  is  usually  defined  as  6,  termed  the  generalized  or  effective 
stepsize,  and  hn  is  the  historical  vector. 

Using  (3.1-6)  to  substitute  for  the  u  in  (3.1-2),  vn  is  found  in 
terms  of  un  and  the  historical  vector 
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8vn  =(8B  +  AM)un  -  AMhn  (3.1-8) 

This  equation  can  be  used  along  with  forms  of  (3.1-6)  for  ti  and  v  to 
remove  everything  except  un,  the  historical  vectors  hn  and  bn,  and 
the  current  forcing  term  Rn  from  (3.1-3).  Making  the  appropriate 
substitutions  leaves 

[m  r 8C  +  82K.]un  =  [m  +  8(C  -  A_1B)]hS  +  8A_1hj;  +  82Rn  (3.1-9) 

Note  3:  Until  now  nothing  has  been  said  about  discretization  by 
means  of  the  popular  second  order  methods,  such  as  Newmark, 
Houbolt,  or  Wilson-0  integration.  These  second  order,  linear 
multistep  methods  have  been  examined  in  detail  by  Geradin  [3.4] 
and  are  easily  cast  in  forms  similar  to  the  first  order  methods.  Two 
integration  approximations  are  required,  one  like  the  first  order  (3.1- 
4),  the  other  in  terms  of  accelerations 

m  m 

X  “iUn-i  =  h  X  Mn-i 

i=0  i=0 

m  .  m 

X  TiUn-i  =  ^  X  Pi“n-i  (3.1-10) 

i=0  i=0 

or  in  simpler  notation 

Un  =  5uun  +  hn 

un  =  hn8uun  +  b£  (3.1-11) 

where  8U  is  8  from  before,  hn  is  as  defined  in  (3.1-7)  and  h„  is 
defined  as 
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K  =  h2Tlbn  -  Cl  =  h2Tl[un-l 

iin-J 

'  Pi  1 

-[“n-1 

M 

•••Un-nJ  1 

1  Pm , 

|  Ym 

(3.1-12) 


Using  (3.1-11)  to  substitute  for  li  and  ii  in  (3.1-1)  and  defining 
5^  =  hri  results  in  a  final  form  similar  to  (3.1-9) 

[m  +  5tiC  +  SiiSuKjun  =  Chn  +  Mhii  +  S^uRn  (3.1-13) 

Since  this  result  is  merely  an  extension  of  the  first  order  form  (by 
choosing  v  =  u  and  using  separate  integration  formulas  for  each),  it 
may  be  treated  in  the  same  fashion  and  is  not  be  dealt  with 
independently. 

3.2  Effects  of  Computer  Implementation 

There  are  two  primary  decisions  to  be  made  when 
implementing  a  numerical  integration  scheme.  The  first  involves 
choosing  one  of  the  equations  already  presented  to  compute  the 
necessary  terms  for  hn.  while  the  second  deals  with  the  choice  of 
weighting  matrices  A  and  B  to  determine  the  auxiliary  vector  v. 
Both  decisions  affect  the  solution’s  stability,  error  propagation,  and 
efficiency. 

3.2.1  Computing  the  Historical  Vector 

The  calculation  of  the  historical  vector  hn  is  referred  to  as  the 
computational  path.  The  three  steps  in  finding  hn,  for  each  of  the 
three  paths  (0,1,2),  are  shown  in  Table  3.2-1.  The  remaining  steps 
which  find  hn,  form  the  right  hand  side  of  (3.1-9),  and  then  solve  for 
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un  are  the  same  for  all  paths  and  are  listed  in  Table  3.2-2.  It  should 
be  mentioned  that  there  are  other  possible  ways  in  which  hn  may  be 
calculated,  but  these  are  the  most  frequently  encountered.  However, 
one  common  variation  is  to  update  iin  within  the  0  path,  creating  a  0’ 
path. 

3.2.2  Choice  of  Auxiliary  Vector 

There  are  two  widely  accepted  choices  for  v.  The  physically 
intuitive  choice  is  to  let  v  =  u.  This  is  called  the  conventional  form  of 
v.  The  other  major  form  was  introduced  in  [3.1]  by  Jensen.  This  form 
is  determined  by 

v  =  Mii  +  Cu  (3.2-1) 

so  that  v=R-Ku.  The  conventional  form  and  Jensen’s  form  are 
numerically  equivalent,  however,  the  effort  and  efficiency  of  each  is 
significantly  different. 

The  conventional  form  has  the  benefit  of  reducing  four  state 
vectors  to  three  ([  u,u,v,v  j  to  [  u,ii,ii  ])  as  well  as  the  physical 
significance  of  v.  Unfortunately,  this  form  sometimes  requires  a 
non-singular  mass  matrix.  Jensen’s  form  does  not  have  this  difficulty 
and  is  more  computationally  efficient  than  the  conventional  form 
(see  [3.2]).  It  should  be  noted,  however,  that  in  a  simulation  coupling 
context  the  auxiliary  vector  and  the  computational  path  are 
determined  by  the  choice  of  integration  package  and  not  selected  for 
efficiency  or  effectiveness  of  implementation  of  the  overall  coupled 
problem. 
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Table  3.2-1  Steps  in  Finding  hn 


Variable 

Equation  Used 

Form  of  Equation 

0  Path 

Vn-1 

(3.1-3) 

Vn-1  +(AC  -  B)tin-1  +  AKun.i  =  ARn_i 

Vn-1 

(3.1-6) 

Vn-1  =  hpoVn-l  +  hn-1 

_ hn _ 

(3.1-7) _ 

,  v  ,.v  v 

_ hn  —  hbn  -  &n _ 

1  Path 

Vn-1 

(3.1-3) 

vn. i  +  (AC  -  B)un_i  +  AKun.i  =  ARn.i 

Vn-1 

(3.1-2) 

Vn-i  =  A  Mikn-i  +  Bun-i 

_ hi _ 

(3.1-7) _ 

_ hn  —  hbn  -  8n _ 

2  Path 

Vn-1 

(3.1-2) 

Vn-1  =  A  Mlin-1  +  Bun-1 

Vn-1 

(3.1-6) 

Vn-1  =(vn-l  -  hn-i)/hPo 

fan 

(3.1-7) 

V  V  v 

hn  =  hbn  -  an 

Table  3.2-2  Remaining  Steps  in  Finding  un 


Variable 

Equation  Used 

Form  of  Equation 

hn 

(3.1-7) 

fcu  -  Uk"  u 

hn  -  hbn  "  Un 

bn 

(3.1-9) 

[m  +  8(C  -  A^B^hS  +  8A_1hn  +  52R„ 

An 

(3.1-9) 

[m  +  5C  +  82k] 

«n 

Solve 

AnUn  =  bn 

_ 

<31-6> _ 

lin  =(un  -  hn)/hPo 
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3.3  Operational  Notation 

Operational  notation  is  presented  to  produce  a  concise 
presentable  notation  form  and  to  facilitate  stability  and  accuracy 
analysis.  The  discrete  Laplace  transform  and  the  z-transform  are 
introduced  to  the  integration  approximation  and  the  historical  terms 
to  produce  an  operational  expression  for  bn,  the  right  hand  side  of 
(3.1-9).  All  of  these  expressions  are  necessary  for  the  evaluation  of 
the  simulation  coupling  algorithm. 

3.3.1  The  Shift  Operator  and  the  Z  Transform 

In  any  series  a  single  term  may  be  related  to  the  previous  or 
following  term  by  means  of  the  shift  operator 

wk  =  M  wk_!  ,wk  =  B'1  wk+1  (3.3-1) 

Repeated  application  will  relate  any  term  wk  to  the  initial  or  final 
term 

Wk  =  Bk  w0  ,  wk  =  Bk'a  wn  (3.3-2) 

Since  the  integration  methods  discussed  so  far  require  only  the 
Wn-i.-.WD  terms,  we  can  focus  on  the  second  of  each  of  the  forms 
above.  Applying  the  discrete  Laplace  transform  to  the  shift  operator 
produces 

wk  =  e(k-n>sh  wn  (3.3-3) 

Substituting  the  standard  z-transform  definition,  z  =  esh,  into  the 
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above  expression 

wk  =  zk‘nwn  (3.3-4) 

With  these  definitions  a  function  f 

f(k)  =  wn  +  aiWn.j  +  a2wn_2  +  a3wn_3  (3.3-5) 

is  transformed  into  the  Laplace  domain 

f(s)  =  (1  +  aie‘sh  +  a2e'2sh  +  a3e'3sh)wn  (3.3-6) 

and  into  the  z  domain 

f(z)  =  (1  +  aiz'1  +  a2z‘2  +  a3z'3)  wn  (3.3-7) 

or  sometimes 

f(z)  =  z‘3(z3  +  aiz2  +  a2z!  +  a3)  wn  (3.3-8) 

Although  for  the  function  listed  above  the  Laplace  and  z  variables 
are  exact,  the  substitutions  which  following  are  only  approximately 
true.  Since  the  time  series  is  obtained  from  numerical  integration 
the  properties  of  the  shift  operator  in  (3.3-1)  and  (3.3-2)  are  not  met 
exactly  as  they  would  be  if  the  true  solution  series  was  used. 

3.3.2  Notation  for  Approximations  and  Historical  Vectors 

Since  numerical  integration  produces  a  time  series  of  vectors 
the  terms  may  be  readily  carried  into  the  Laplace  or  z  domains. 
Given  an  integration  approximation  of  the  form  previously 
introduced  (where  ao  is  normalized  to  1  and  now  factoring  (3o  out, 
leaving  p'i  =  pi/po) 
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m  m 

w n  +  Z  aiWn-i  =  8  wn  +  2  P'iWn-i 


"  n  1  ^  v*1’  n-i  w  n  1  ^  n-ij  (3  3  9) 

Applying  the  shift  operator  and  taking  the  discrete  Laplace 
transform  of  the  approximation  gives 


or  in  the  z  domain 


X  aie*iShjwn  =  fi|l  +5)  P,ie'iShjwn 

+  X  “iz "‘jwn  =  5^1  +  X  P'iz"‘jw,i 


(3.3-10) 


(3.3-11) 


By  rearranging  (3.1-6)  the  historical  vector  may  be  written  as 


h  n  —  Wn  -  5\V  n 


(3.3-12) 


Using  this  definition  for  hn,  (3.3-10)  and  (3.3-11)  may  be 
manipulated  to  form  the  following  expressions: 

C(S)  =  s(£  ^e-i5jw„-(£aie-ish)ffn 
h"(z)  =  a|X  P'iz ’‘jwn  -  |X«iz’‘|wn  (3.3-13) 


It  is  useful  to  define  commonly  occurring  polynomials  in  z  and  esh 


ui 

P(.)  =1  +  2  ®i(») 


a(.)=l  +  lP'i(.) 


(3.3-14) 


Then  in  either  domain 
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p(.)wn  =  8o(.)wn  (3.3-15) 

This  also  allows  the  historical  vector  to  be  written  as 

hn(.)  =  8(a(.)-l)wn  +  (l-p(.))  wn  (3.3-16) 

Another  useful  operator  definition  which  is  needed  is 

wn  =[p(.)/8a(.)]wn  =  \y(.)  wn  (3.3-17) 

3.3.3  Notation  for  Forcing  Term  bn 

The  presence  of  hn  in  the  forcing  term  requires  that 
computational  path  be  dealt  with  in  order  to  come  up  with  the 
correct  notation.  Although  each  path  has  a  distinct  operational  form, 
the  choice  of  auxiliary  vector  has  no  effect.  Recalling  the  basic  form 
of  bn 

bn  =  [m  +8(C-A'1B)]hn  +  8A‘1hn  +  82Rn  (3.3-18) 

.  •  .  u 

Using  the  definitions  of  the  previous  section  hn  can  be  written 

bn(.)  =  (l-8\|/(.))un  (3.3-18) 

V 

which  is  also  independent  of  path.  The  form  of  hn  for  each  path 
comes  from  the  equations  in  Table  3.2-1  and  when  combined  with 
the  above  equations  the  following  general  form  results: 

bn  =  +  <|>cC+  <|>kK]  un  +  <j>R.Rn  (3.3-19) 

These  constants  are  detailed  in  Table  3.3-1. 
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Table  3.3-1  Operational  Constants  in  bn 


Path 

<|»M 

tyc 

0K 

0 

1  -5\p 

0 

82-8v1 

Sy1 

i 

1  -(p2/a) 

S(l-p) 

S2(l-a) 

s2 

8  a 

2 

1  -(by)2 

8(l-8¥) 

0 

1 

3.4  Use  Of  Predictors 


Simulation  coupling  is  made  possible  by  the  ability  to 
accurately  predict  wn  from  earlier  knowledge  of  w.  A  general  linear 
multistep  predictor  is  suggested  in  [2.6]  for  use  in  large  coupled 
problems.  The  form  of  this  predictor  is 


m 


m 


Wn  =  X  a^Wn-i  +  hpo  X  PiWn-i  +  (hpo) 

i=l  i=l 


,2  Jj2,  p  . . 

X  YiWn-i 


i=l 


(3.4-1) 


As  with  the  integration  methods  it  is  also  beneficial  to  place  the 
predictor  in  operational  form  by  defining  polynomials  in  esh  and  z  as 
follows: 

p  -i 

Pw(.)  =  X  ai  (•) 
i=l 

Pw(.)  =  s£pf(.)"  (3.4-2) 

i=l 

2  £1,  p  -i 

Pw(.)  =  5  X  Ti  (-) 
i=  1 
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Using  these  definitions  as  well  as  (3.3-16)  gives  the  final  form  of  the 
predictor  as: 

Wn(.)  =  p(.)wn  =  [pw(.)  +  8\J/pw(.)  +  (5\j02pw(.)]wn  (3.4-3) 


3.5  Examples 

3.5.1  Trapezoidal  Rule 

The  trapezoidal  rule  is  a  one  step  (m=l)  implicit  method 
commonly  used  in  many  integration  applications.  Its  constants  are 

a0  =  1,  ai  =  -1,  p0  =  Pi  =  0.5 

Then  the  operational  forms  are 

P(esh)  =  1  -  e"sh  ,  p(z)  =  1  -  z'1 
a(esh)  =  1  +  e'sh  ,  cr(z)  =  1  +  z"1 

y(esh)  =2.1..-  e~ —  -2.  tanh(sh/2),  \|/(z)  =  2. =  2. 

h  !  +  e-*h  h  h  x  +  z-i  h  z  +  1 

and  the  historical  forms  are 

hn(s)  =  (0.5  hwn+  wn)e"sh 

hn(z)  =  (0.5  hwn+  wn)z’ 

3.5.2  Gear’s  Two  Step  Method 

Another  implicit  method  introduced  by  Gear  in  [3.5]  is  a  two 
step  A-stable  form.  Its  constants  are 

ao  =  1,  ai  =  -4/3,  CX2  =  1/3 
Po  =  2/3,  P!  =p2  =  0 

The  operational  forms  are 
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p(esh)  =  1  -  Vsh  +  V2s\  p(z)  =  1  -  4z-i+lz-2 
3  3  3  3 

o(esh)  =  1  =  a(z) 

¥(e>h)  =(J_)(l  -  +  = 

\2h/l  3  3/  2  h 

¥z)={M'-fz'+3z)=^2h4yl 

and  the  historical  forms  are 

hn(s)  =  ^(4e"Sh  "  e’2Sh)wn 

hn(z)  =  ^4z_1  -  z’2)wn 
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Chapter  Four 


Constraint  Equations  and 
Numerical  Simulation  Coupling 


Before  the  simulation  coupling  process  may  be  fully  developed 
and  analyzed,  the  incorporation  of  constraints  into  equations  of 
motion  must  be  covered.  The  remaining  portion  of  this  chapter 
provides  the  details  of  simulation  coupling.  Finally,  the  analysis  of 
the  stability  and  accuracy  of  the  coupling  process  is  addressed. 

4.1  Methods  of  Handling  Constraints 

In  an  attempt  to  develop  an  accurate  and  reliable  method  of 
dealing  with  constraints,  many  different  procedures  have  been 
suggested  [4. 1-4.6].  Two  methods  of  primary  interest  for  use  in 
coupled-field  problem  are  penalty  formulations  and  more  traditional 
Lagrange  multiplier  forms.  Besides  the  basic  implementations  of  the 
methods  there  are  additional  stabilization  procedures  to  improve 
performance  which  are  covered.  Both  methods  detailed  here  are 
developed  from  using  traditional  variational  principles  and 
specifically  the  Euler-Lagrange  equations. 
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4.1.1  The  Lagrange  Multiplier  Method 

A  typical  n-DOF  mechanical  system  is  defined  by  the  following 
energy  products: 

Kinetic:  T  =  ^-tiTMii  and  Potential:  V  =  4TKu  (4.1-1) 

2  2 

as  well  as 

r=  Cii  (4.1-2) 

where  T  are  Rayleigh's  dissipative  forces,  representing  damping 
proportional  to  velocity.  Additionally,  the  response  of  the  system  is 
restricted  by 

d>(u,ii,t)  =  0  (4.1-3) 

where  d>  is  a  vector  of  p  individual  constraints. 

When  no  derivative  terms  are  present,  these  functions  may  be 
used  to  eliminate  p  DOFs  in  the  vector  u.  Unfortunately,  this  is  can  be 
difficult  and  it  is  not  always  clear  which  individual  degrees  should  be 
removed.  Instead,  Lagrange's  “method  of  the  undetermined 
multiplier”  may  be  used  to  form  a  term  which  represents  the  forces 
which  ensure  the  constraint  is  satisfied.  This  term,  along  with  the 
nonconservative  externally  applied  forces,  8W  =  R8u,  and  energy 
products,  may  be  combined  to  form  the  action  integral  for  a 
constrained  system 

Jr‘2 

(l  +  W  +  <t>TX  +  rT  u)  dt  (4.1-4) 
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where  X  is  the  vector  of  Lagrange  multipliers  and  L  is  the  Lagrangian 
T  -  V.  By  applying  the  principle  of  least  action,  8 A  =  0,  the  Lagrange 
equation  of  motion  is 


d_  W  .  .  r  -  ~  X  =  R 

Ldt\auj  du  /  3u 


(4.1-5) 


By  using  (4. 1-1, 2)  in  (4.1-5)  the  matrix  equation  of  motion  for  a 
constrained  system  is  obtained 

Mii  +  Cii  +  Ku  =  R  +  GtX  (4.1-6) 

where  G  is  the  Jacobian  of  (4.1-3)  with  respect  to  u. 

By  taking  the  derivative  of  twice  forms  the  following 
equation: 

Gii  +  Gti  =  0  (4.1-7) 

Using  (4.1-6)  to  substitute  for  the  second  derivative  in  the  above 
equation  leaves  a  form  which  is  solvable  for  X. 

GM  1GTX=  GM ^[Cii  +  Ku  -  R]  -  Gii  (4.18) 

Since  M  is  large,  taking  the  inverse  is  a  costly  operation  To  correct 
this  the  state  vector  is  partitioned  into  the  degrees  of  freedom  which 
contribute  to  the  constraints  and  those  which  do  not,  u  =  [ufT  ucT]T. 
Generally  uc  is  small  compared  to  Uf.  If  a  lumped  mass  model  is 
used  then  the  following  substitutions  may  be  made  in  (4.1-8): 

Gt  =  gIc  M_1=M^  G  =  Glc  (4.1-9) 
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where 


GT=[  o  GtJ  M'*  = 


Mf/  0 

0  Mil 


G  = 


0  G, 


CC  j 


With  these  definitions  the  Lagrange  multipliers  are 
Xnh  =  [GccMcUjJc]  GccMCc[Cii  +  Ku  '  R]  -  [GccMciGic]  Gccii  (4.1-10) 


where  the  inverses  are  now  taken  of  matrices  with  the  dimension  of 
ucc  and  dimension  p. 

holonomically  constrained  problems  may  be  treated  in  much 
the  same  way.  Taking  the  first  derivative  of  <E>  when  derivative 
terms  are  present  leaves 

Hti  +  Gti  =  0  (4.1-11) 

where 

H  =  ^ 

dll 

This  equation  is  the  same  in  structure  as  (4.1-7)  with  H  and  vj 
substituted  for  G  and  G.  Making  the  same  assumptions  as  the 
nonholonomically  constrained  problem,  the  Lagrange  multipliers  for 
the  holonomic  problem  are 

Kh  =  [hccMccH Jc]" 1 H CCM cc [C ii  +  Ku  -  R]  -  [hccMc1chJ'c]'1GcCu(4.1-12) 

There  is  one  problem  with  this  formulation.  Using  (4.1-7)  and 
(4.1-11)  forces  the  second  and  first  derivatives  of  to  be  zero  and 
not  <1>  itself.  Because  of  this  fact  the  values  of  X  need  to  be  corrected 
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by  some  form  of  stabilization  technique.  At  this  time  a  single  best 
way  to  do  this  does  not  exist  (see  [4.5]). 


4.1.2  The  Penalty  Method 

As  derived  in  the  previous  section,  the  matrix  equation  of 
motion  for  a  holonomically  constrained  system  is 

Mu  +  Cii  +  Ku  =  R  +  GTA 


The  above  equations  have  n+p  unknowns  in  u  and  A..  The  penalty 
method  provides  the  additional  p  equations  by  defining  A  in  the 
following  manner: 


A  =  -k<&,  as  - — »  0 

K 


(4.1-13) 


In  general  the  penalty  constant  k  does  not  need  to  be  the  same  in 
each  of  the  p  constraints,  in  which  case  k  becomes  a  p  x  p  diagonal 
matrix.  Substitution  into  (4.1-6)  leaves 

Mii  +  Cii  +  Ku  =  R  -  Gtk<&  (4.1-14) 

The  penalty  method  format  implicitly  assumes  that  the  constraints  of 
(4.1-3)  are  violated. 

One  difficulty  in  this  type  of  formulation  should  be  noted. 
Although  this  method  converges  to  the  correct  solution,  once  an  error 
occurs  there  is  no  way  in  which  the  energy  associated  with  the 
penalty  correction  may  be  dissipated.  To  correct  for  this  excess 
energy,  stabilization  procedures  can  be  introduced.  For  the  penalty 
method  one  such  stabilization  routine  is  detailed  here  (see  also  [4.5]). 
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It  is  useful  to  note  that  the  penalty  formulation  of  (4.1-14)  can 
be  achieved  by  augmenting  the  Lagrangian  with  a  fictitious  potential 
representing  the  compromised  constraints 


V*  =  d>T  a  k  $ 
2 


(4.1-15) 


To  dissipate  energy,  fictitious  Rayleigh's  forces  are  included 


f  =-ad& 
dt 


(4.1-16) 


T  .  *  * 

Replacing  the  term  d>  A,  with  T  -V  in  the  action  integral  and  again 
applying  the  principle  of  least  action 


WaL).|L)^=R.^TJdd>+K<, 

\dt\du  3ui  du  Ut 


(4.1-17) 


and  making  the  appropriate  substitutions 


Mii  +  Cii  +  Ku  =  R  -  GTce(^  +  k0> 

\dt 


(4.1-18) 


where  both  a  and  k  may  be  p  x  p  diagonal  matrices.  In  this  form  a 
acts  as  the  penalty  and  k  becomes  the  decay  constant  for  the  error  in 
the  constraints. 

Although  the  terms  introduced  above  apply  to  the  holonomic 
case,  it  is  possible  to  derive  the  system  equations  of  motion  for  the 
nonholonomic  case  in  a  similar  manner.  When  derivatives  are 
present  in  the  constraints  a  fictitious  kinetic  energy  and  a  fictitious 
set  of  Rayleigh’s  dissipative  forces  similar  to  (4.1-15)  and  (4.1-16) 
augment  the  Lagrangian  as  follows: 
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T*  =  ^$TaK$ 

2 


(4.1-19) 


and 


T  =  -  a  d> 


(4.1-20) 


Placing  these  terms  within  the  action  integral  yields: 


MaL).  9L  _r  =R.a^1a^ 
4t\3u/  ^u;  3u  Idt 


(4.1-21) 


Making  the  appropriate  substitutions  to  obtain  this  equation  in 
matrix  form  results  in 


Mu  +  Cii  +  Ku  =  R  -  HT 


(4.1-22) 


4.2  The  Simulation  Coupling  Process 

The  solution  of  dynamic  problems  by  means  of  simulation 
coupling  requires  three  distinct  elements.  The  first  of  these  elements 
is,  obviously,  the  individual  simulations  for  the  separate  fields. 
Additionally,  elements  to  evaluate  and  correct  for  any  violation  of 
the  constraints  which  couple  the  fields  and  to  handle  data 
management  tasks  and  control  the  execution  of  the  separate 
simulations  are  necessary. 

There  are  two  different  ways  in  which  the  execution  of  the 
coupling  process  may  take  place.  The  simplest  form  starts  with  a 
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Prediction  -  Corrective  Force 
Evaluation 
I  n  =  n  +  1  T”™ 


Interaction  Variables 


Figure  4.2-1  Information  Flow  in  Concurrent  Evaluation 

prediction  of  the  field  states  (x  and  y)  in  order  to  produce  a  set  of 
corrective  forces  based  on  the  incorporation  of  constraints.  Then  the 
individual  simulations  calculate  the  field  responses  at  the  given  time. 
Finally,  time  is  incremented  and  any  necessary  updates  to  the  data 
structure  are  handled.  This  sequence  is  illustrated  in  Figure  4.2-1.  In 
this  case  the  evaluations  made  by  the  simulations  may  be  carried  out 
in  series  or  in  parallel.  This  execution  order  is  called  concurrent 
evaluation. 

The  second  execution  process  is  slightly  more  involved  than  the 
concurrent  evaluation.  It  also  starts  with  a  prediction  of  the  state 
variables  and  a  corrective  force  prediction.  However,  only  one  field 

evaluation  is  carried  out.  This  execution  is  used  to  update  the  force 

evaluation  before  analyzing  the  second  field.  Again,  the  process  ends 
with  time  update  and  data  management  functions.  This  is  called  a 
staggered  evaluation  and  is  shown  as  Figure  4.2-2.  In  staggered 

evaluation  one  field  may  be  consistently  analyzed  first  or  the  order 

may  be  switched  as  seems  fit. 
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Figure  4.2-2  Information  Flow  in  Staggered  Evaluation 

As  introduced  in  Chapter  Two,  the  dynamics  of  a  two  field 
coupled  system  can  be  expressed  as 

Mxx  +  Cxx  +  Kxx  =  Rx  +  fcx(x,x,y  ,y) 

Myjr  +  Cyy  +  Kyy  =  Ry  +  fcy(x,x,y,y)  (4.2-1) 

where  the  corrective  forces  due  to  the  constraints,  fcx  and  fcy,  are 
found  using  the  methods  already  discussed  and  other  similar 
methods.  For  the  following  discussion  the  forces  are  assumed  to  be 
calculated  using  a  stabilized  penalty  formulation.  It  is  useful  to  note 
that  these  equations  are  uncoupled  if  fcx  and  fcy  are  considered  to  be 
unresolved  externally  applied  forces. 

Solving  these  equations  by  means  of  the  numerical  integration 
was  discussed  in  Chapter  Three,  the  system  equations  become 

[mx  +  6XCX  +  SjKjjxn  =  [<|>MxMx  +  <|>cxCx  +  (^KxKJ  Xn  +  $Rx(Rx,n  +  fcx,n)(4 .2-2) 
[My  +  SyCy  +  8yICyj  y  n  =  [<|>MyMy  +  <(»CyCy  +  <|>KyKy]yn  +  <|*Ry  (R  y,n  +  ^cy,n) 

or  in  simpler  notation 

[(l-4>M)Mu  +  (8-<j)c)Cvl  +  (82-(|>K)Ku]un  =  <|>R(Rn  +  fc,n)  (4.2-3) 

where 
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Mx 

0 

Cu  = 

cx 

0 

r 

n 

Kx 

0 

Rn  = 

Rx,n 

0 

My 

9 

0 

Cy  . 

> 

0 

Ky. 

9 

.  Ry,n_ 

S  -[S^Syf  _<>,=[  f  _  fC,0  =  [fCX>„fCy,n]r  _  1  =  [1  If 

Although  the  forms  of  the  constraining  force  developed  in 
Section  4.1  do  not  require  linear  constraints  be  used,  analysis  of  the 
simulation  coupling  process  is  greatly  aided  by  a  linear 
approximation  as  follows: 


=  Ccii  +  Kcti  (4.2-4) 

With  such  an  approximation  the  forces  of  constraint  are 

Holonomic  ^c>n  =  'Kc<xKc(iin  +  Kun)  (4.2-5) 

Nonholonomic  fc?n  =  -Cja(Cciin  +  KcUn) 

The  manner  in  which  these  forces  are  evaluated  is  based  entirely  on 
the  particular  execution  process  used. 

4.2.1  Concurrent  Evaluation  Coupling 

As  has  already  been  mentioned  concurrent  evaluation  applies 
the  same  correction  term  to  each  field.  For  this  term  to  be  treated  as 
an  applied  force  (allowing  the  simulations  to  remain  uncoupled),  it 
cannot  be  a  function  of  the  current  time  step.  Through  use  of  the 
predictors  introduced  in  Chapter  Three,  any  dependence  on  the 
current  states  is  removed. 

First  the  integration  approximation  is  applied  to  remove 
derivative  terms 
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Holonomic 
Nonholonomic  fc,n 


(4.2-6) 


or  in  operational  form 


Holonomic  fc>n  =  -Kj  aKc(y(.)I  +  k)ud  (4.2-7) 

Nonholonomic  fc>n  =  -Cj  a(Cc\|f(.)2  +  KcV(.))un 


In  this  form  the  base  states,  un,  are  estimated  using  a  predictor  of 
the  form  of  (3.4-1).  The  predicted  correction  terms,  after 
substituting  p(.)un  for  u£,  are 

Holonomic  fc.n  =  -KjaKc(\|/(.)I  +  K)p(.)un  (4.2-8) 

Nonholonomic  f£n  =  -Cj  a(Ccy(.)2  +  Kc¥())p()un 

Making  the  following  definitions  to  simplify  the  notation 

M^h  =  (v(-)I  +  k)p(.)  (4.2-9) 

M-nh  =  (l¥(-)2  +  C'c1Kc\|/(.))p(.) 

Then  the  system  equation  for  the  concurrent  evaluation  is 

[(1  -<1>m)Mu  +  (8-<j>c)Cu  +  (82-<J>k)Ku  +  4>r  KjaKcM-h] “n  =  <!>r  Rn  (4.2-10) 
or  for  the  nonholonomic  case 

[(1  -4m)Mu  +  (8-<j>c)Cu  +  (82-<^k)Ku  +  <(>r CjaCcpJ un  =  $r Rn  (4.2-11) 


The  implementation  of  the  concurrent  evaluation  is  given  as 
Table  4.2-1. 
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Table  4.2-1  Implementation  of  Concurrent  Simulation  Coupling 


Step  in  Coupling  Process  Associated  Equation  Number 


1.  Use  states  n-1  to  n-m  to  find  u£  (3.4-1) 

2.  Use  states  n-1  to  n-m  to  find  hn  and  hn  (if  needed)  (3.1-7) 

3.  Use  u£,  hn,  and  hn  to  find  fp  n  (4.2-6) 

4.  Send  fpn  to  separate  simulations 

5.  Separate  simulations  solve  for  un  (4.2-3) 

6.  Calculate  tin  and  iin  if  simulations  do  not  provide  them  (3.1-6) 

7.  Update  time  and  increment  n  =  n+1 _ 


4.2.2  Staggered  Evaluation  Coupling 

The  staggered  evaluation  follows  directly  from  the  concurrent 
evaluation.  However,  since  the  fields  must  be  handled  in  series,  the 
response  from  the  first  field  is  used  instead  of  a  prediction  when 
calculating  the  correction  force  for  the  second.  The  operational  form 
also  comes  from  the  concurrent  evaluation  in  the  following  manner. 


fp  - 

*c,n  — 


form  of  the  correction  force,  fc,n. 

(kJoKcL  (KlaKck, 

(KjaKcU  (kJaKcix 

H  Jn 

(4.2-12) 


or 


fp 

‘c,n 


fealdi  (kJokJ 
0 


*y 


0 


{kJoKcV*  (kJoikJ 

The  predictor  term’  is  expanded  as  follows: 


yx 


K] 

i 

i 

.  yn  . 

(4.2-13) 
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(4.2-14) 


The  upper  half  of  this  matrix  deals  with  the  prediction  of  each  field 
for  evaluating  the  constraint  forces  acting  on  the  first  field;  the  lower 
half  does  the  same  for  the  second  field.  By  eliminating  the  term  in 
the  lower  half  which  corresponds  to  the  prediction  of  the  first  field, 
the  actual  values  are  used  instead  of  predicted  ones.  This  results  in 
the  following  matrix: 

Ph  0 

0  Ph 

(y(.)I  +  ic)  0  (4.2-15) 

0  Ph 
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The  implementation  of  the  staggered  evaluation  is  given  in 
Table  4.2-2. 

Table  4.2-2  Implementation  of  Staggered  Simulation  Coupling 


Step  in  Coupling  Process  Associated  Equation  Number 


1.  Use  states  n-1  to  n-m  to  find  u£  (3.4-1) 

2.  Use  states  n-1  to  n-m  to  find  hn  and  (if  needed)  (3.1-7) 

3.  Use  u£,  and  h„  to  find  fc,n  for  x  field  (4.2-6) 

4.  Send  fc,n  to  x  field  simulation 

5.  Simulation  solves  for  xn  (4.2-3) 

6.  Use  yn,  h£,  hjj,  and  xn  to  find  fc,n  for  y  field  (4.2-6) 

7.  Send  fc,n  to  y  field  simulation 

8.  Simulation  solves  for  yn  (4.2-3) 

9.  Calculate  tin  and  iin  if  simulations  do  not  provide  them  (3.1-6) 


4.3  Analysis  of  the  Simulation  Coupling  Process 

4.3.1  Stability  Analysis  of  Simulation  Coupling 

Having  developed  operational  expressions  for  the  system 
equation  of  motion,  the  stability  analysis  is  fairly  straight-forward. 
There  are  two  possible  methods  to  accomplish  the  analysis.  Both 
consider  the  unforced  response  (Rn  =  0  for  (4.2-11))  and  develop  a 
system  characteristic  matrix.  Since  the  response  is  unforced  the 
eigenvalues  of  the  matrix  must  be  less  than  or  equal  to  unity,  to  keep 
the  response  from  growing  with  time.  Although  both  methods  will 
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be  shown  in  detail  in  the  examples  of  Chapter  Five,  they  will  be 
briefly  introduced  here. 

The  first  method  uses  the  matrix,  C  (z),  which  is  the  matrix 
multiplying  un  in  (4.2-11).  For  instance,  the  characteristic  matrix  for 
the  concurrent  evaluation  and  nonholonomic  constraints  is 

C(z)  =  (l-<t>M)Mu  +  (8-<J)c)Cu  +  (82-(j>K)Ku  +  <j>RKcaKc|xni1  (4.3-I) 

Taking  the  determinant  of  this  matrix  produces  the  characteristic 
polynomial,  C(z).  The  simulation  coupling  process  is  stable  if  the 
roots  of  C(z)  lie  within  the  unit  circle.  Optimally,  it  is  desired  that 
these  roots  reflect  the  same  stability  as  the  fully  coupled  system.  An 
example  of  this  would  the  trapezoidal  rule;  an  optimally  stable 
coupling  process  would  keep  this  method’s  A-stable  characteristics. 

The  second  method  relies  on  a  time  domain  analysis  and  is  a 
common  numerical  analysis  procedure  (see  [4.7]).  The  state  variable 
un  can  be  written  as  a  function  past  state  variables  and  derivatives 
using  a  form  of  (3.1-9) 

[m  +  8C  +  82k]uii  =[m  +8(C-A'1B)]h}j  +  8A'1hn  +  82fc,n  (4.3-2) 

Using  the  integration  approximation  for  u  and  (4.3-2)  to  substitute 
for  un,  the  derivative  un  can  also  be  written  as  function  of  past 
states.  The  definition  of  the  auxiliary  vector  and  integration 
approximation  for  v  provide  vn  and  v  in  terms  of  the  past  states. 
These  equations  make  it  possible  to  write  a  matrix  which  when 
multiplied  by  the  vector  [  un.i  vn_i  un-2  vn-2  •••  Un-m  vn.m  ]  will  give 
the  vector  [  un  vn  un-i  vn-i  ...  un.m+i  vn.m+i  ].  The  roots  of  this 
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vector  must  also  be  less  or  equal  to  one  for  the  coupled  simulation  to 
be  stable.  Although  developing  this  matrix  is  more  difficult  than 
finding  C(z),  then  eigenvalues  may  be  calculated  directly  from  this 
point  without  first  taking  the  determinant  to  find  the  characteristic 
polynomial. 

Unfortunately  using  these  methods  is  not  nearly  as  easy  as  it 
would  seem.  Large  coupled  problems  tend  to  have  thousands  of 
degrees  of  freedom,  and  the  number  of  roots  to  the  characteristic 
equation  can  be  several  times  that  number.  This  makes  calculation 
of  these  roots  difficult  even  with  computer  methods.  Additionally, 
the  only  parameters  the  dynamicist  has  available  to  ensure  stability 
are  the  integration  time  step,  the  choice  of  the  predictor,  and  the 
penalty  used  to  find  the  corrective  force.  The  choice  of  overall 
penalty  is  generally  governed  by  requirements  to  keep  the 
constraints  violations  very  small  and  the  time  step  must  be  small 
enough  to  stabilize  the  high  frequency  poles  associated  with  the 
application  of  a  penalty.  This  leaves  only  the  predictor  and  at  this 
time  there  are  no  guaranteed  choices  for  optimally  stable  predictors. 

4.3.2  Accuracy  Analysis  of  Simulation  Coupling 

The  accuracy  analysis  also  involves  manipulating  the 
characteristic  equation,  but  this  time  the  s  domain  form  is  desired. 
By  expanding  the  exponential  terms  in  the  operational  expressions, 
e*sh  =  1  -  sh  +  0.5(sh)2  -  ...  ,  the  characteristic  equation  may  be 
written  in  powers  of  s.  This  allows  the  characteristic  equation  in  s  to 
be  transferred  back  to  a  differential  expression  and  from  there  to  a 
standard  eigenvalue  problem.  As  an  example  consider  an  undamped 
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system  integrated  using  the  trapezoidal  rule  and  computational  path 
2.  The  characteristic  equation  is 


[(fiyfM  +  52k]ud  =  0 


(4.3-3) 


Substituting  for  the  operational  expressions  from  Chapter  Three 


1  +  e_sh 


M  +  K  un  =  0 


(4.3-4) 


Expanding  e*sh 


Jl  -  ±{sh)2  +  0(sh)4jM  s2  +  K  un  =  0 


(4.3-5) 


This  will  give  the  differential  equation 


Mii  --Lh2M'ii'  +0(u<4>)+Ku  =0 
6 


or  an  eigenvalue  problem 


+  lh2K  +  ...)(<d h)2  -  O(coh)4  +  K  =  0 


(4.3-6) 


(4.3-7) 


These  equations  show  that,  as  expected,  the  trapezoidal  rule 
introduces  no  artificial  damping  but  does  produce  a  certain  amount 
of  frequency  distortion.  The  same  sort  of  analysis  may  be  carried 
out  on  the  simulation  coupling  characteristic  equation.  The 
numerical  values  can  then  be  compared  to  similar  values  obtained 
for  the  fully  coupled  system. 
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Chapter  Five 


Illustrative  Numerical  Examples 


Although  the  sii  .  .lation  coupling  theory  is  aimed  at  the 
solution  of  large  scale  problems,  it  is  not  possible  to  show  all  the 
details  of  analysis  of  such  problems.  For  this  reason  the  following 
chapter  is  organized  to  present  a  series  of  smaller  dimensioned 
problems  in  which  the  all  of  the  aspects  of  the  theory  presented  are 

shown.  The  first  example  is  the  one  dimensional  motion  of  a  single 

rigid  body.  The  response  is  considered  subject  to  different  forcing 

conditions.  The  second  example  is  a  more  complicated  model 

representing  a  satellite  held  at  the  end  of  the  shuttle’s  manipulator 
arm.  This  example  shows  some  of  the  effects  of  imposing  constraints 
on  a  flexible  domain.  Finally,  the  single  rigid  body  is  again 
considered.  All  six  degrees  of  freedom  are  allowed  to  show  the 
interaction  between  translational  and  rotational  motion  in  the 
presence  of  constraints.  All  examples  use  the  concurrent  evaluation 
along  with  a  stabilized  penalty  format  to  find  the  constraint  forces. 
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5.1  One  Dimensional,  Single  Body  Example 

As  a  single  body  example  consider  the  structure  illustrated  in 
Figure  5.1-1.  This  body  is  representative  of  a  dual  spin  stabilized 
satellite  (see  [5.1]).  The  mass  properties  for  each  section  are 
Rotor  =  700  kg  Platform  =  300  kg 

Each  center  of  mass  (CM)  is  located  on  the  spin  axis.  The  platform 
CM  is  0.5  m  up  from  the  joint  between  the  bodies  and  the  rotor  CM  is 
0.75  m  down  from  the  joint.  The  complete  center  of  mass  is  0.375  m 
down  from  the  joint. 

The  inertia  properties  of  the  satellite  about  their  CMs  are: 

Spin  Transverse 

Rotor  350.0  kg-m2  306.25  kg-m2 

Platform  150.0  kg-m2  100.0  kg-m2 

Complete  500.0  kg-m2  734.375  kg-m2 
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Figure  5.1-1  Satellite  Rigid  Body  Model 
The  satellite  is  subject  to  the  following  forcing  conditions 
resulting  from  attitude  control  type  maneuvers: 

Forces  =  1  N  any  axis  Torques  =  10  N-m  about  spin 

The  forces  and  torques  are  applied  on  the  rotor  section.  For  the 
purpose  of  simulation  coupling  the  single  body  is  divided  into  two 
separate  rigid  bodies  where  the  platform  joins  the  rotor. 


5.1.1  Analysis  of  1-D  Example 

Considering  the  response  of  each  body  in  the  direction  of  the 
spin  axis,  the  equations  of  motion  are 

Mrxr  =  Rr  and  Mpxp  =  Rp  (5.1-1) 


Additionally,  the  bodies  are  rigidly  connected  along  the  spin  axis,  or 
xr  -  xp  =  0.  By  applying  the  stabilized  penalty  method  discussed  in 
Chapter  Four,  the  following  constraint  forces  are  found: 


fc=-[l  -l]Ta[(xr-Xp)  +  K(xr-xp)] 
Adding  these  forces  to  the  system  leaves 

'  Mr  0  ljkrj  /  Rr 

.  0  Mp  J\xp/  |  Rp 


(5.1-2) 


(5.1-3) 


Subtracting  the  second  equation  from  the  first  and  rearranging 
terms  produces  the  constraint  error  equation 


Mint  v  \  _ _ _M 

ip  Mrivip 

The  penalty  k  controls  the  steady  state  error 


(xr-xp)  -  a-^-(xr-Xp)  +  aK-^-(xrxp; 

JVirMp  MrMr 


Rj..  Rp. 

Mr  Mp  (5.1-4) 


aK^-|(xr-xp)sJ 


MrM 


RrMp  -  RpM, 


MrMp 


(5.1-5) 
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With  the  given  mass  properties  and  a  requirement  that  the  steady 
state  error  be  less  than  1.0x10* 7  meters 


k=J--3.0x105  .  _  , 

a  (5.1-6) 

The  other  penalty,  a,  determines  the  damping  ratio  of  the  error.  For 

a  ratio  of  0.3,  a  =  2400.0  which  leaves  k  =  125.0. 

The  equations  of  motion  may  be  integrated  following  the  (0) 

computational  path  and  using  the  trapezoidal  rule 

Wn-Wn.!  =0.5h(wn  +  wn.i)  (5.1-7) 

Using  this  equation  to  substitute  for  the  accelerations  in  (5.1-3)  gives 


xr 

Xp 


n-1 


Mr 

0 


0 

M 


P  -I 


Rr 


R 


+  (fc)n 


p  ;n 


(5.1-8) 


The  constraint  force  is  found  using  a  last  value  predictor,  xpn  =  xn_!, 
leaving 


fc,n  — 


1 

-1  ‘ 

■ 

Xr\  .  1 Xr 1 

=  -a 

.  -1 

1  . 

'  +  K{  r 

Up /n-1  ^P/11-1. 

(5.1-9) 


These  two  equations  are  used  to  carry  out  the  simulation. 

To  perform  the  analysis  of  the  simulation  coupling  the  notation 
introduced  in  Chapter  Three  is  used.  Defining  u  =  [  xr  xp  ]T,  Rr  =  0, 


Rp  =  0,  and  M  = 


Mr 

0 


0 

M 


leaves 


un  =  un.!  +  hun-i  +-^— un-i  -  cc^-M 

4  4 


-ll 


1  ‘l][Un-I  + 
-11, 


KUn-l] 


(5.1-10) 
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Using  the  characteristic  equation  for  holonomically  constrained 
coupled  problems  listed  in  (4.2-10) 

[(Hm)M  +  o*rKJkc(¥  +  K)p(z)]u„  =  0  (5.1-11) 

For  this  example  the  last  value  predictor  is  characterized  by  the  shift 
operator,  z*1,  and  the  polynomials  <)>M  and  4>r  are  listed  in  Table  3.3-1. 
Making  the  appropriate  substitutions  results  in  the  following  form: 

[(8y)M  +  a^xir^KjKcCy  +  kK1]^  =  0  (5.1-12) 


The  operational  notation  for  the  trapezoidal  rule  is  found  in  Section 
3.5.  Using  the  given  polynomials  and  clearing  most  fractions  leaves 


z(z-1)2M  +  ^kJkJ(z-1)(z+1)  +  £k(z+l)2 
2  L  2 


Un  =  0 


(5.1-13) 


From  the  form  of  the  constraint,  <&  =  xr  -  xp,  the  matrix  Kc  is  found  to 
be  [1  -1]T.  Expanding  the  polynomials  and  matrices  gives  the 
characteristic  matrix  equation,  C(z) 


C(z)  = 


(z3-2z2+z) 


Mr  0  +  tth-f  1 "  ^[(l+^b-jz2  +  (Kh)z  +  f—  - 

0  MD  J  2  L-1  1  JL\  2  /  '2 


(5.1-14) 


Defining  a  temporary  variable,  a(z)  =  (l+^-jz2  +  (ich)z  +  (^L-l 


C(z)  = 


Mr(z3-2z2+z)  +  c^-a(z) 
-  fi^-a(z) 


-  ah-a(z) 

2 

Mp(z3-2z2+z)  +  c^-a(z) 


(5.1-15) 


For  stability,  the  roots  of  the  determinant,  det(C),  must  be  within  the 
unit  circle.  Taking  the  determinant 
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z(z-l)2  z3  +  </Ml2Ll^  +  ffiKh.j  .  2jz2  +  |/MX<x  jgKh  Ll}2  +  //M2oL|gK_h 

UMeMpA  2  4  1)  UM,MpA  2  /  f  4  2  /j 

(5.1-16) 

Since  the  first  three  roots  (z  =  0,1,1)  do  not  vary  with  a,  k,  or  h, 
they  are  not  as  important  as  the  remaining  three.  For  the  penalties 
given.  Figure  5.1-2  shows  two  of  the  remaining  three  roots  moving 
away  from  z=l  as  the  step  size  is  varied  from  0  to  0.008,  where 
0.008  is  at  the  point  of  instability.  The  third  pole  moves  from  the 
origin  toward  a  value  of  z  =  0.02.  The  desired  time  step  should  be  as 
large  as  possible  while  still  being  stable.  For  this  reason  h  is  chosen 
to  be  0.005. 


Figure  5.1-2  Pole  Locations  of  Rigid  Body  Model 

The  choice  of  a  does  not  seem  to  have  a  great  effect  on  the 
overall  stability.  If  a  is  doubled  the  limiting  time  step  increases  by 
less  than  5%.  Therefore,  since  a  is  proportional  to  the  damping  ratio 
and  large  damping  ratios  mean  small  overshoot,  it  is  often  feasible  to 
increase  a  by  as  much  as  a  factor  of  five.  Additionally,  changing  the 
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mass  ratio  of  platform  to  rotor  does  not  have  any  effect  on  the 
coupled  simulation  stability.  On  the  other  hand  k  has  a  great  effect 
on  stability.  Doubling  the  value  of  k  will  almost  double  the  stable 
time  step. 

5.1.2  Response  of  1-D  Example 

The  response  is  considered  subject  to  three  separate  forcing 
conditions;  a  constant  force  input,  a  ramp  input,  and  a  pulse  sequence 
input.  In  each  case  the  responses  of  the  rotor  and  platform  are 
compared  to  the  response  of  the  complete  rigid  case  for  a  two  second 
run.  The  time  step  used  is  h  =  0.005  and  the  penalties  are  those 
suggested  in  the  previous  section.  The  position,  velocity,  and 
acceleration  for  the  constant  forcing  case  are  illustrated  in  Figure 
5.1-3. 

The  error  between  the  coupled  and  complete  response  is 
noticeable  only  in  the  acceleration  plot  and  then  only  during  the  first 
two  seconds.  The  steady  state  error  settles  out  to  exactly  1.0  x  10* 7 
as  desired.  To  check  the  calculated  stability  boundary,  a  coupled 
simulation  was  made  with  h  =  0.008.  The  constraint  error  and  the 
constraint  error  rate  are  shown  in  Figure  5.1-4. 

The  second  forcing  condition  is  a  ramp  thrust  from  0  to  1  N 
over  a  two  second  period.  The  position,  velocity,  and  acceleration  for 
this  case  are  presented  in  Figure  5.1-5.  In  this  case  the  constraint 
error  also  increases  with  time  reaching  a  final  value  of  1.0  x  10‘7. 
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Time 

Figure  5.1-3  1-D  Example  Response  to  Constant  Forcing 
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The  third  forcing  condition  is  a  series  of  1  N  thrust  pulses.  Each 
pulse  lasts  a  quarter  second  and  separated  from  the  next  pulse  by 
the  same  interval.  The  position,  velocity,  and  acceleration  for  this 
case  are  illustrated  in  Figure  5.1-6.  Some  noticeable  differences 
occur  in  this  case  due  to  the  rapidly  changing  input.  However,  the 
coupled  position  and  velocity  still  track  the  actual  solution  fairly 
well. 
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Time 


Time 

Figure  5.1-5  1-D  Example  Response  to  Ramped  Forcing 


Figure  5.1-6  1-D  Example  Response  to  Pulsed  Forcing 
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5.2  One  Dimensional  Shuttle-Satellite  Example 

The  next  example  is  the  one  dimensional,  four  mass  model  of  a 
satellite  held  in  the  space  shuttle  manipulator  arm  illustrated  in 
Figure  5.2-1.  In  actuality  the  connection  between  the  end  effector  of 
the  arm  and  the  satellite  platform  is  assumed  rigid,  and  this  will 
provide  the  coupling  between  the  separate  fields  of  the  shuttle 
dynamics  and  the  satellite  dynamics.  The  shuttle  partition  is  made 
up  of  the  shuttle  and  the  manipulator  arm.  The  one  dimensional 
idealization  of  the  flexibility  in  the  arm  is  a  spring  connecting  the 

two  masses.  The  satellite  partition  is  similar  to  the  satellite  in  the 

previous  rigid  example  with  the  addition  of  a  small  amount  of 
flexibility  between  the  rotor  and  the  platform.  The  masses  involved 
are 

Shuttle  Body  (Ms)  85,000  kg 

Manipulator  Arm  (Ma)  140  kg 

Satellite  Platform  (Mp)  300  kg 

Satellite  Rotor  (Mr)  700  kg 

The  flexibilities  involved  are 

Manipulator  Flexibility  (Ka)  300  N/m 

Satellite  Flexibility  (Krp)  1800  N/m 


Shuttle  Partition  Satellite  Partition 


Rigid  Connection 
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Figure  5.2-1  1-D  Shuttle-Satellite  Model 

Forces  are  allowed  to  the  system  from  the  control  jets  on  the  shuttle 
and  the  satellite.  The  satellite  jets  are  capable  of  the  same  1  N  thrust 
while  the  shuttles  jets  can  deliver  up  to  2  kN  (see  [5.2]  and  [5.3]). 


5.2.1  Stability  Analysis  of  Shuttle-Satellite  Example 
The  equation  of  motion  for  each  partition  is 


‘  Ms 

0 

(ks)  + 

-Ka 

(X«) 

0 

Ma  . 

Wal 

L-Ka 

Ka  J 

Ual 

and 


Krp  -Krp  jxp  |  /  0  \ 

-Krp  Krp  iXri  \Rr| 


(5.2-1) 


The  complete  system  equation  of  motion  found  by  eliminating  a 
degree  of  freedom  due  to  the  rigid  connection,  xa  -  xp  =  0 


Ms  0  0 

0  Ma+Mp  0 
0  0  Mr 


lxS| 

;xp 

IXri 


Ka 

Ka 

0 


-Ka 


0 


Ka+Krp  -Krp 


-K 


rp 


K 


rp 


' 

fxs 

Rs 

< 

Xp 

0 

Urj 

Rr 

(5.2-2) 


The  equivalent  simulation  coupled  system  with  corrective  constraint 
forces  is 


Ms  0 
0  M, 
0  0 
0  0 


0 

0 

Xs  \ 

0 

0 

X.  \ 

+ 

Mp 

0 

Xp 

0 

Mr 

Xr  | 

Ka  -Ka  0  0 

-K,  Ka  0  0 

0  0  Krp  -Krp 

0  0  -Krp  Krp 


where  the  constraint  force  is 


X$ 

Xa 

RS  \ 
fc  1 

Xp 

( 

-fc  j 

Xr 

l  Rr  I 

(5.2-3) 


fc  =  -a[(xa-xp)  +  K(xa-Xp)] 


8  1 


(5.2-4) 
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To  determine  the  penalty  constants,  a  and  k,  an  error  analysis 
similar  to  the  one  used  in  the  first  example  is  done.  Subtracting  the 
third  line  of  (5.2-3)  from  the  second 


Ma  +  Mp  ^ 
MaMp  C 


(5.2-5) 


Rearranging  terms,  defining  e  =  xa  -  xp,  and  substituting  for  fc 


e  +  a 


Ma  +  M£ .  +aKM1±MEe  =  Kffi(X[  _  .  ik(xa .  x$) 


M,Mr 


MaMp  Mp 


Ms 


(5.2-6) 


Replacing  the  terms  on  the  right  hand  side  with  terms  from  the  first 
and  fourth  lines  of  (5.2-3)  leaves 

”  +  gMa  +  Mp  ^  +  aKMa  +  Mp  ^  _  (Rr  -  Mrxr)  _  (Rs  -  Msxs) 

MaMp  MaMp  Mp  Ma  (5.2-7) 


Using  (5.2-7)  to  require  the  steady  state  error  be  less  than  10' 6 
meters 


K  a1’510*  (5.2-8) 

The  other  penalty,  a,  determines  the  damping  ratio  of  the  error.  For 
a  ratio  of  0.1,  a  =  24,000.0  which  leaves  k  =  6,250.0. 

Since  using  C  (z)  to  do  the  stability  analysis  would  require 
taking  the  determinant  of  a  4x4  matrix  of  third  order  polynomials, 
the  second  method  of  determining  stability  is  used  instead.  Applying 
the  trapezoidal  rule  and  a  last  value  predictor,  (3.1-9)  becomes 

[m  +  82K]un  =  82Mun-i  +  25Miin-i  +  Mun.j  +  82fc,n  (5.2-9) 
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where  M  and  K  are  the  defined  in  the  coupled  equation  (5.2-3). 
After  substituting  for  fc>n 

M  +  82K]un  =  52Miin-i  +[25M  -a52G]iin.i  +[m  -a52GK]un.!  (5.2-10) 

where  G  =  k£Kc  =  [  0  1  -1  0]T[  0  1  -1  0].  Using  the  the  trapezoidal 
rule  iin  may  be  written 

tin  =  6_1(Un  -  un_!)  -  lin-1  (5.2-11) 

or  using  (5.2-10)  for  un 

M  +  82K]un  =  8Mii„.i  +  [m  -82K  -a8G]iin-i  +[-5K-a8GK]un.!  (5.2-12) 
Doing  the  same  thing  for  iin 


M  +  82K]iin  =  [-82K]iin.i  +  [-28K  -  aG]un-i  +  [  -K-  aGic]un.  j  (5.2-13) 

Defining  KEff=|M  +  82Kj  and  using  (5.2-10,-12,-13)  leaves  the 
following  amplification  matrix  relating  states  at  tn-i  to  states  at  tn: 


ke;I-62k_ 

Keh[-28K  - 

ccG] 

K^-K-ccGk] 

KfJ,8M 

-  62K 

-a8G] 

KeII  -8K-  oSGk] 

(5.2-14) 

Ke!,«2M 

K^l28M  -  a82Gj 

¥^\[m  -  q82Gk] 

All  of  the  eigenvalues  of  this  matrix  must  be  within  the  unit  circle  to 
ensure  the  coupling  process  is  stable.  Performing  this  analysis,  the 
maximum  time  step  that  maintains  stability  is  h  =  1.50x1  O'4. 

5.2.2  Accuracy  Analysis  of  Shuttle-Satellite  Example 

The  characteristic  equation  for  an  undamped  structural  system 
with  holonomic  constraints  is 
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Hm)M  +  (82-<)>k)K  +  cc<|>rG(y  +  ic)p(sh)Jun  =  0  (5.2-15) 

For  the  0  computational  path  using  the  trapezoidal  rule  and  a  last 
value  predictor  the  following  constants  are  defined: 

4>m  =  l-8v>  <(>k  =  8  -SxjT1,  <|)r  =  Sxj/"1,  p(sh)  =  esh,  y  =  tanhj^j 

With  these  constants  (5.2-15)  becomes 

Stanhj^-jM  +  8tanh  +  aStanh  ^^GOanh^  +  K)esh  un  =  0  ^  ^ 

or  after  clearing  fractions 

tanh^^-jM  +  K  +  aG(tanh/s^Lj  +  K)esh  un  =  0 

Using  the  series  expansion  for  tanh  and  e 

tanh(^)=  s(l  -^-s2h2+— s4h4-  ]  and  esh=  (l+sh+J-s2h2+i-s3h3+^-s4h4+  I 
12/  \  12  240  I  \  2  6  24  / 

which  makes  the  complex  terms 

tanh 2(SJ2-)  =  s2[l  -^-s2h2+-U-s4h4-  •) 

V  2  /  \  6  720  / 

tanhfe-kje511  =  s(l+sh+-5-s2h2+~/-s3h3+-i—  s\4+  ) 

V  2  /  \  12  12  240  / 

Inserting  these  terms  into  (5.2-17)  leaves 
Ms2(l  -lsV+Ji-sV--)  +  K 

\  6  720  / 

+  aGs(l+sh+^-s2h2+-Ls3h3+ 

V  12  12 

+  aG  Kll+sh+-/-s2h2+-/-s 
\  2  6 

or  rearranging  terms 
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|M--^-Ms2+aGh+aG^-s2+0(h4)  )s2 
+  (aG(l+Kh)  +aG(^-h2+  ^-)s2+  0(h4))s 
+  (K-HxGK+^aGKs2+0(h4)) 


un  =  0 


(5.2-19) 


Dropping  higher  order  terms  and  using  the  fact  that  [Ms2  +  K]un  =  0 


M-^-K-HxGh+aG^-M^K  s2 
\  6  12  / 

+  jaG(l+Kh)+aG(^-h2+ k|-.)M‘1k)s 
+  |k+oGk+  hi  aG  kM  1 K 


U  n  =  () 


(5.2-20) 


This  equation  can  be  used  to  find  the  coupled  system  frequencies. 
Table  5.2-1  compares  the  continuous  system  frequencies  to  the 
frequencies  shifted  by  integration  alone  (discrete  frequencies)  and 
the  frequencies  shifted  by  integration  and  coupling.  As  (5.2-20) 
shows  some  artificial  damping  occurs  as  a  result  of  the  simulation 
coupling  process.  However,  the  damping  ratios  of  the  coupled 
frequencies  is  less  than  10-7.  The  is  also  a  coupled  frequency 
corresponding  to  the  constraint  penalties  at  approximately  1200 
rad/sec. 


Table  5.2-1  Comparison  of  Frequency  Shi 

ting  Effects  (rad/sec) 

Continuous 

Discrete 

Coupled 

0.00000000 

0.00000000 

0.00000000 

0.50034156 

0.50034155 

0.50034119 

2.66408400 

2.66408398 

2.66408308 
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5.2.3  Response  of  Shuttle-Satellite  Example 

Two  different  responses  were  considered  for  this  example. 
Both  cases  involved  forcing  the  system  using  the  shuttles  control  jet 
thrusters.  Each  run  is  compared  to  the  identical  conditions  using  the 
complete  system  shown  in  (5.2-2).  For  both,  the  complete  and  the 
coupled  systems,  a  time  step  of  0.0001  seconds  was  used.  The 
penalties  suggested  in  Section  5.2.1  were  used  for  the  coupled  case. 

The  first  case  is  a  20  second  run  where  the  force  of  2  kN  has 
been  applied  for  the  first  second.  Figure  5.2-2  shows  the  position  of 
the  shuttle,  the  satellite  rotor,  and  the  joint  between  the  arm  and  the 
platform.  Additionally,  Figure  5.2-3  shows  the  constraint  error  and 
constraint  error  rate  for  this  case. 

Since  this  is  an  undamped  example  the  system  energy  should 
be  constant.  However,  due  to  the  stabilization  of  the  constraints, 
energy  will  be  bled  off  from  the  system.  This  process  will  continues 
as  long  as  there  is  any  potential  energy  stored  due  to  flexibility.  This 
artificial  damping  is  a  very  slow  process  and  has  little  noticeable 
effect,  but  care  should  still  be  taken  to  make  the  damping  ratio  for 
constraints  less  than  or  equal  to  the  damping  present  in  the 
structure,  if  possible. 

The  second  case  should  show  that  no  significant  energy  is  lost 
even  over  a  period  of  many  time  constants.  It  is  a  two  minute  run 
where  the  2  kN  force  is  applied  for  ten  seconds.  The  figures  for  this 
run  shows  a  twenty  second  period  at  the  end  of  the  two  minutes. 
Figure  5.2-4  shows  the  translations  of  the  major  system  parts. 
Figure  5.2-5  illustrates  the  error  between  the  coupled  response  and 
the  complete  response  for  the  shuttle,  the  rotor  and  the  end  effector. 
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Figure  5.2-2  Shuttle-Satellite  Example  Short  Term  Response 
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Complete 
Coupled  -  Arm 
Coupled  -  Platform 


Figure  5.2-4  Shuttle-Satellite  Example  Long  Term  Response 


5.3  Free  Floating  Rigid  Body  Example 

The  final  simple  example  is  a  return  to  the  satellite  described 
in  Section  5.1  and  illustrated  in  Figure  5.1-1.  This  example,  however, 
will  consider  all  six  DOFs,  R0  =  [  X  Y  Z  ]T  and  0  =  [  8i  02  03  F-  The 
translations  express  displacement  from  the  inertial  frame  in  inertial 
coordinates  and  the  angles  represent  Euler  angles  of  a  2-1-3  rotation 
from  the  body  frame  to  the  inertial  frame.  The  equations  of  motion 
for  a  rigid  system  with  these  specified  DOFs  are  derived  in  many 
dynamics  texts  (see  [5.4]). 

For  the  fully  rigid  case  there  are  six  active  constraints.  Three 
translational  constraints  require  the  point  where  the  rotor  joins  the 
platform  to  have  the  same  inertial  position.  Three  rotational 
constraints  require  the  Euler  angles  for  the  platform  orientation  and 
the  rotor  orientation  to  be  equal.  The  two  sets  of  constraints  are 
distinct  enough  to  justify  having  a  unique  a  and  k  for  the  translations 
and  the  rotations. 

5.3.1  Analysis  of  Rigid  Body  Example:  Translational  DOFs 

Each  body  frame  is  fixed  at  the  body  CM  with  the  z  axis 
oriented  along  the  spin  axis  indicated  in  Figure  5.1-1.  If  ra  is  the 
position  of  the  joint  between  the  platform  and  the  rotor  with  respect 
to  the  CM  of  frame  a  in  frame  a  coordinates.  With  this  definition  the 
constraints  on  the  translational  DOFs  are 

<I>  =  Ro  -  R£  +  jCb  rr  -  !Cg  rP  (5.3-1) 
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where  !Cb  is  the  matrix  representing  an  orthogonal  2-1-3  rotation 
from  body  to  inertial  frames.  Taking  the  first  two  derivatives  gives 

^  =  Rr„  -  RS  +  'fir'  -  >CV  +  'Cefo/xr')  -  >c|(o>PxrP)  (J  3  2) 

and 

^f-  =  Ro  -  RS  +  rcirr  -  +  ^(co'xr1)  -  ^(©PxrP) 

dt  (5.3-3) 

+  ICB(cbrxrr)  -  rCg(cbpxrP)  +  ^^(©^(ci^xr1-))  -  TCg(coPx(ci)PxrP)) 


These  terms  are  used  to  form  the  constraint  correction  force,  fc,  and 
the  constraint  error  equation. 

The  translational  equations  of  motion  for  these  two  bodies  are 


R«=^,CBRr+„bf‘  a"d  ro=S7F 


(5.3-4) 


Taking  the  difference  of  these  equations  and  substituting  e  for  d>  and 


-oc[^  +  k<I> 
L  dt 


for  fc  leaves 


Ro  -  R?  =  —  iCgRr  -  —  a[e  +  ice] 
u  °  mr  a  rarmp  1  J 


(5.3-5) 


Adding  the  additional  terms  to  form  e  and  the  constraint  error 
equation 


e  +  m:-+3«c  +  a  Ke  =  ±K?BRT  +  f* 


m  rlQ  p 


HI  rlQp 


mr 


(5.3-6) 


where  f*  includes  the  nonlinear  terms  treated  as  applied  forces 

f*  =TCBrr  -  !C^rP  +  IC^((Drxrr)  -  IC«(o)PxrP)  + ICB(cbrxrr) 

-  !C |(cbpxrP)  +  ICB(corx(corxrr)) -  ^^(©Px^PxrP))  (5.3-7) 
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In  order  to  keep  the  analysis  linear,  the  nonlinear  forcing  term 
is  assumed  to  be  much  smaller  than  the  actual  forcing.  This  leaves 
the  vector  form  of  the  error  equation  found  in  (5.1-4).  Given  the 
same  error  limits  and  desired  damping  (steady  state  error  less  than 
10*7  and  damping  of  30%),  the  penalty  constants  and  time  step  limit 
found  in  the  first  example  will  still  hold.  These  numbers  are 

a  =  2,400.0  k  =  125.0  h  <  0.008  (5  3.8') 


5.3.2  Analysis  of  Rigid  Body  Example:  Rotational  DOFs 

As  described  earlier  the  rotational  constraints  require  the  Euler 
2-1-3  angles  for  each  frame  to  be  equal 

<I>  =  0r  -  6P  (5.3-9) 


Taking  the  first  two  derivatives  gives 


d$.  =  er  .  ep  =  s;1©1  -  SpV 


dt 


(5.3-10) 


and 


d!^  =  e'r .  e'p  =  s-r1d)r-  Sp^  +  s'rV-  Sp©p  (5.3-1 1) 


-1  .  --i 


dt 


where  Sr  is  a  non-orthogonal  transformation  from  Euler  angle  rates 
to  angular  velocity  coordinates  (see  [5.4]).  These  terms  are  used  to 
form  the  constraint  correction  torques,  xc,  and  the  constraint  error 
equation 

The  rotational  equations  of  motion  for  these  two  bodies  are 
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d)r  =  If  1 T  r  -  Ir  1(d)rxIrCDr)  +  I^lS 
(bp  =  -Ip  1(copxIpcop)  -  Ip'1-!? 


(5.3-12) 


Premultiplying  by  the  appropriate  form  of  Sa1  and  taking  the 
difference  of  these  equations,  as  well  as  substituting  e  =  and 


= -Saa 


do  +  ko 

dt 


leaves 


S'rV  -  SpV  =  Sr’1(Sr+  SpjS^I-r^Ir  +  IpjfpVPxIpW15-^1^) 

-  S'r1(Sr  +  Sp)Sp1Ir1(lr+  IpJlp^Sr+SpJafe  +  Ke]  +  Sr‘ 1  I^T  r(5 -3‘ 1 3 > 


Adding  the  additional  terms  to  form  e  for  the  constraint  error 
equation  and  making  the  following  definitions:  S  =  Sr ^Sr  +  SpjSp1  and 
I  =  Ir^Ir  +  Ipjlp1 ,  leaves 

e  +  S I(Sr  +  Spjae  -f  S I{Sr  +  Sp)a ice  =  S^V  +  T*  (5.3-14) 
where  X*  includes  the  nonlinear  terms 

X*  =  S l(coPxIpO)P  -  CDrxIrtor)  +  s/co1  -  Sp^P  (5.3-15) 

Two  assumptions  are  made  to  keep  the  analysis  linear.  First,  the 
nonlinear  forcing  term,  X*,  is  assumed  to  be  small  with  respect  to  the 
actual  forcing.  Additionally,  since  the  error  between  the  Euler  angles 
is  forced  to  be  small,  the  transformation  matrix  for  the  two  body 
frames  is  assumed  to  be  the  same,  Srp.  This  will  reduce  the  matrix  S 
to  2S‘rp.  Rearranging  the  remaining  linear  terms  leaves 

S'rJj^Srpe  +  4  a  e  +  4cuce  =  S^Ip(lr  +  Ip)'1! r  (5.3-16) 
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Using  this  equation  and  desired  values  for  steady  state  error  and 
damping  ratio,  the  penalty  constants  can  be  found  in  the  same 
manner  as  the  first  and  second  examples.  With  a  steady  state  error 
of  10‘4  degrees  and  and  damping  ratio  of  0.1  the  constants  are  a  = 
375.0  and  k  =  40.0. 

To  find  the  maximum  stable  time  step  the  complete  equations 
of  motion  must  be  formed.  Combining  (5.3-12)  and  the  fact  the  Sa 
relates  Euler  angle  rates  to  angular  velocity,  the  complete  first  order 
unforced  equation  of  motion  can  be  formed 

0  0 
Srp  0 

0  0 
0  0 
0  It1 
0  0 
0  -Ip1 
0  0 

or  redefining  terms 

0  =  K0  -  X*  -  Ga[©  +  k©]  (5.3-18) 


0 

0 

0 

0 


Ir‘ 


0 

Ip1 

0 


0 

0 

tor 

0 

0 

er 

1 

0 

0 

©p 

| 

Srp 

0 

,  ep  , 

\ 

S  rptt 


I'rWxlKO1) 

0 

Ip  *(  ft)PxIpC0P) 
0 


•  r  \ 

(0  1 

(0r  \ 

•  r  | 

6  1 

er  1 

i 

}  +  K1 

) 

“P 

0)P 

ep  ) 

0P  )_ 

(5.3-17) 


Using  the  trapezoidal  rule  and  a  last  value  predictor,  an  amplification 
matrix  similar  to  the  one  in  (5.2-14)  can  be  formed.  With  Keff 
defined  as  I(12)-8K,  the  matrix  is 


K'e^SK  -  aG]  KeIJk  -  ccGk] 
KE1fl8I(12)  -  8aG]  KE1Jl(12)  -  8ccGk] 


(5.3-18) 
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By  requiring  the  eigenvalues  of  this  amplification  matrix  to  be  within 
the  unit  circle  the  critical  time  step  is  found  to  be  h  =  0.03. 

5.3.3  Response  of  Rigid  Body  Example 

For  the  response  tne  penalties  used  are 

a  k 

Translation  2400.0  125.0 

Rotation  375.0  40.0 

The  maximum  time  step  is  0.008  seconds  and  the  actual  time  step 
used  is  0.005  seconds.  The  satellite  (rotor  and  platform)  is  given  an 
initial  spin  rate  of  0.5  rpm.  The  complete  set  of  plots,  all  degrees  of 
freedom  and  errors,  is  given  in  Appendix  A. 

The  first  case  reflects  a  constant  force  input  of  1  N  in  the  z  axis 
Due  to  the  satellite’s  symmetry  this  forcing  only  excites  the  z  axis 
translation.  The  error  response  is  exactly  as  shown  in  the  first 
example,  with  constraint  error  damping  out  at  two  seconds  to  the 
steady  state  value  of  10"7  meters.  The  second  case  is  similar  with  a 
constant  torque  input  of  10  N-m  about  the  z  axis.  Again  only  the  63 
DOF  shows  any  response.  This  response  damps  out  after  four  seconds 
to  a  steady  state  value  of  10'4  degrees.  These  responses  show  that 
the  simulation  will  perform  as  desired  under  the  simplest  of  forcing 
conditions. 

Although  the  earlier  cases  verify  performance  of  the  coupling 
algorithm  the  more  interesting  cases  arc  those  designed  to  show  the 
cross-coupling  effects  between  the  translation  and  rotation  DOFs. 
The  third  case,  a  1  N  force  along  each  axis  shows  one  new  effect.  The 
constant  spin  about  the  third  axis  produces  a  long  period 
displacement  error  in  the  x  and  y  axes.  This  displacement  is  a 
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maximum  at  1.5xlO'6  meters,  larger  than  the  desired  value  but  still 
small  enough  that  the  penalty  constants  do  not  need  to  be  changed. 
Had  this  value  not  been  sufficient  the  constants  may  be  changed 
easily.  For  instance,  if  the  maximum  allowable  error  is  10*6  meters 
then  by  raising  k  by  a  factor  of  1.5  and  a  by  Vl .5  then  maximum  is 
reduced  while  still  still  maintaining  reasonable  damping. 
Additionally,  if  the  penalties  are  changed  by  a  significant  amount  it 
may  also  be  necessary  to  adjust  the  time  step. 

The  final  case  involves  applying  a  0.1  second  1  N-m  torque  to 
one  of  the  non-spinning  axes.  In  this  case  the  nonlinear  terms  which 
were  ignored  in  the  analysis  of  Section  5.3.2  start  to  have  a 
significant  impact  on  the  motion.  The  original  constants  must  be 
altered  to  account  for  this  nutation  of  the  spinning  satellite  and 
ensure  the  stability  of  the  simulation  coupling  process.  Running 
several  different  penalty  combinations,  simply  adding  a  large 
amount  of  damping  seems  the  most  effective  correction  for  nutation 
effects.  After  achieving  stability,  the  penalties  are  modified  to  keep 
the  steady  state  errors  small.  The  total  of  these  two  corrections 
raises  a  by  a  factor  of  forty  and  k  by  a  factor  of  four.  The  time  step 
is  dropped  during  this  process  to  0.0025,  again  to  ensure  stability. 

There  is  one  final  effect  of  a  rigid  body  simulation  to  discuss. 
When  using  Euler  angles  it  is  possible  to  encounter  singularities 
where  the  three  angles  condense  to  one  DOF.  This  continues  to  be  the 
case  in  a  simulation  coupling  process.  The  most  successful  way  to 
correct  for  this  effect  is  to  switch  to  a  different  Euler  angle 
combination  when  close  to  a  singularity.  Although  this  will 
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successfully  avoid  the  singularity  problems,  for  this  example  the 
initial  conditions  were  altered  instead  to  avoid  any  singularities. 

5.4  Conclusions  from  Simple  Examples 

These  examples  should  show  that  for  small  cases  it  is  possible 
to  write  an  error  equation  to  find  appropriate  penalties  and  to 
complete  a  linear  stability  analysis  to  find  a  time  step  for  the 
simulation  coupling  algorithm.  Both  rigid  and  flexible  problems 
display  unique  difficulties  to  watch  for  as  larger,  more  involved 
problems  are  attempted. 

For  rigid  cases  the  most  difficult  problem  to  handle  is  the 
singularities  Each  simulation  handles  them  differently  and  it  may 
be  necessary  to  change  constraints  and  even  penalties  to  account  for 
the  singularity.  The  only  other  noteworthy  problem  is  how  to  deal 
with  the  coupling  between  displacement  and  rotation  DOFs  which 
may  occur  due  to  rotating  frames,  nutation,  or  simple  geometry 
considerations.  However,  by  making  small  adjustments  in  the 
penalty  constants  and  time  step  these  difficulties  can  usually  be 
handled. 

Simulation  coupling  of  flexible  cases  presents  the  unique 
problem  of  an  additional  dissipative  force  in  the  penalty  stabilization 
scheme.  Care  should  be  taken  not  to  significantly  alter  the  system 
damping  ratio,  otherwise  the  constraint  error  will  damp  out  but  the 
response  will  not  converge  to  the  correct  solution. 

Finally,  both  flexibility  and  rotation  effects  produce  significant 
long  period  error  responses.  This  long  period  motion  is  usually  much 
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larger  than  predicted  steady  state  maximum  values  and  a 
should  always  be  made  for  this  sort  of  constraint  error  mode. 


check 
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Chapter  Six 


Space  Station  Example 


The  Space  Station  Assembly  Complete  model  used  in  this 
example  consists  of  3852  degrees  of  freedom  and  1577  elements.  It 
is  made  up  of  three  primary  bodies:  the  center  span  and  two 
articulating  “wings”.  The  individual  bodies  are  separated  at  the 
alpha  joints  pictured  in  Figure  6.1  and  the  gamma  and  beta  joints  are 
modelled  as  fixed.  Since  it  is  possible  to  set  up  the  complete  three- 
body  multibody  problem  for  comparison  the  Space  Station  makes  a 
good  test  case  for  the  validity  of  simulation  coupling  as  an 
alternative  simulation  strategy. 

For  the  purpose  of  simulation  coupling  the  response  of  each 
body  is  calculated  according  to  the  single  flexible  body  equations 
presented  in  [6.1].  One  individual  rigid-flex  one-body  simulation  is 
used  for  each  of  the  separate  bodies.  The  simulation  makes  use  of  a 
popular  fourth  order  Runge-Kutta  integration  scheme  instead  of  the 
linear  multi-step  methods  which  are  suggested  in  Chapter  Three. 
Additionally,  the  Space  Station  is  too  large  for  a  flexible  analysis  like 
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Port  PV Arrays  Starboard  PV  Arrays 

8159  8134  8115  8150 


Figure  6.1.  Space  Station  Freedom  Assembly  Complete 


the  one  used  in  Section  5.2  to  be  done.  Because  of  these  two  factors 
the  approach  taken  is  slightly  different  than  the  examples  of  the 
previous  chapter. 

For  the  non-articulated  simulations  all  6  DOFs  at  each  alpha 
joint  are  rigidly  constrained.  The  articulated  simulations  allows  free 
rotation  about  the  y  axis  (02  degree  of  freedom).  The  penalties,  a 
and  k  are  to  be  designed  to  keep  the  translation  constraint  violation 
less  then  10  8  feet  and  the  rotation  constraint  violation  less  than  10- 6 
degrees. 

6.1  Analysis  for  the  Space  Station  Example 

As  has  been  mentioned  the  flexible  model  of  the  Space  Station 
is  too  large  for  a  full  analysis.  It  is  still  reasonable,  however,  to  do 
the  rigid  body  analysis.  This  will  allow  a  first  guess  at  the  complete 
model  penalties  using  the  rigid  body  penalties.  If  this  guess  should 
prove  unacceptable  then  further  adjustments  can  be  made  based  on 
the  error  responses. 
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The  analysis  of  the  Space  Station  is  broken  into  two  two-field 
problems:  the  interaction  between  the  left  wing  and  center  span  and 
the  interaction  between  the  right  wing  and  center  span.  The 
constraint  error  equations  for  a  rigidly  coupled  two-field  problem 
derived  in  Section  5.3  are 

e  +  me  +  mw  a  g  +  mc  +  rnw.a  Ke  _  JLiCrRc  +  f*  ,,  .  . . 

Translation  mcmw  mc  (o-l-l) 

and 

Rotation  e  +  S I(SC  +  Sw)ae  +  S I{SC  +  Sw)a ice  =  Sc^Ic1*0  +  x*  (6.2-2) 

were  center  (c)  and  wing  (w)  properties  and  attendant  matrices 
replace  the  previous  rotor  and  platform  ones. 

Using  these  equations  and  again  assuming  the  nonlinear  terms 
are  much  smaller  than  the  applied  forcing  the  penalty  constants  for 
both  the  starboard  and  port  wings  are 

a  k 

Translation  90,000  8,000 

Rotation  1,150,000  90.0 

The  translational  penalties  limit  the  time  step  at  0.0001  seconds. 

6.2  Response  for  the  Space  Station  Example 

Due  to  the  the  type  and  location  of  the  reaction  jet  clusters 
there  are  two  forcing  conditions  commonly  used  for  comparison  of 
the  Space  Station  model.  They  are  a  one  second  100  lb  force  in  the  z 
direction  and  a  0.2  sec  torque  about  the  x  axis  of  2500  ft-lb.  Using 
the  two  forcing  conditions  there  are  three  cases  to  be  considered:  a 
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rigid  non -articulated  case,  a  rigid  articulated  case,  and  a  flexible  case 
with  articulation.  The  plots  for  all  cases  are  shown  in  Appendix  B. 

Running  the  rigid  non-articulated  case  shows  two  areas  where 
the  penalties  must  be  adjusted  to  ensure  the  proper  response.  First, 
the  penalties  on  the  rotations  must  be  raised  substantially  to  counter 
the  coupling  which  the  inertia  matrix  provides.  Overall,  k  is  raised 
by  a  factor  of  eight  and  a  is  doubled  before  the  error  response  is 
acceptable.  Also,  the  non-primary  responses  need  additional 
damping.  Primary  responses  are  the  DOF  being  directly  forced  and 
the  articulation  DOF,  which  corresponds  to  02-  By  overdamping  the 
other  error  responses  coupling  between  error  terms  is  virtually 
eliminated.  The  overdamping  involved  raising  a  by  additional  factor 
of  ten  and  lowering  k  by  a  factor  of  twenty. 

One  numerical  effect  of  simulation  coupling  seems  to  be  the 
creation  of  alternate  response  paths  which  solve  the  error  equations 
but  do  not  provide  the  correct  solution  for  the  coupled  equations  of 
motion.  For  the  Space  Station  one  of  these  paths  is  the  unconstrained 
response  of  the  center  span.  Under  certain  penalty  values  the  error 
will  decay  rapidly  to  zero,  leaving  the  center  span  ‘blind’  to  the  wing 
on  either  side.  In  fact,  the  basic  penalties  listed  in  Section  6.1 
converged  to  the  unconstrained  solution. 

After  correcting  the  penalties  the  other  rigid  cases  have  little 
new  information  to  offer.  There  are  a  few  small  errors  between  the 
multibody  solution  and  the  simulation  coupling  solution  but  all  are 
within  acceptable  limits.  The  constraint  errors  on  the  average  are 
much  smaller  than  specified  due  to  the  overdamping  effects. 
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The  flexible  model  includes  10  flexible  modes  for  each  single 
body.  The  included  modes  are  the  first  ten  based  on  frequency. 
When  switching  to  the  flexible  case  the  penalties  do  not  need  to  be 
changed  significantly  but  the  time  step  is  cut  in  half  to  account  for 
the  flexibility.  Unfortunately,  the  same  accuracy  seen  in  the  rigid 
cases  is  not  maintained  in  the  flexible  cases. 

The  flexible  case  shows  mosts  of  the  errors  expected  form 
simulation  coupling.  There  are  some  noticable  frequency  shift  effects 
as  well  as  significant  artifical  damping  effects  present.  These  effects 
can  normally  be  reduced  by  tightening  the  penalties  used.  In  the 
simple  example  of  Section  5.3  these  effects  were  not  serious  because 
the  error  was  kept  small.  In  this  case  the  penalties  cannot  be  raised 
to  a  sufficent  level  to  counter  these  effects  for  two  reasons.  First, 
pushing  the  penalties  higher  seriously  limits  the  time  step,  making  it 
unreasonably  small.  Also,  the  penalties  needed  to  remove  most  of 
these  effects  ends  up  causing  overflow  problems  for  the  computer  on 
which  the  code  was  written.  Some  of  these  difficulties  may  be  solved 
using  the  numerical  improvements  discussed  in  Appendix  C. 

The  largest  difference  in  the  flexible  cases  appears  as  a 
frequency  shift  in  the  lowest  frequency  mode.  In  the  x  torque  case, 
for  instance,  this  shift  is  most  evident  in  the  63  DOF.  This  error  also 
causes  a  slow  drift  in  the  other  angular  DOFs  as  well  as  reinforcing 
the  low  frequency  motion  in  the  translational  DOFs.  The  same  sort  of 
error  is  also  present  in  the  z  force  case.  Also,  the  z  force  case 
appears  unstable,  however,  the  actual  response  shows  the  first  mode 
grow  in  magnitude  for  ten  seconds  before  it  begins  to  decay. 
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Chapter  Seven 


Conclusion 


Simulation  coupling  is  a  procedure  which  allows  individual 
existing  simulations  of  separate  fields  to  be  combined  for  the  solution 
of  coupled  field  problems.  It  is  a  potentially  time  and  cost  saving 
method  for  dealing  with  these  problems.  If  the  constraints  which 
couple  the  separate  fields  are  adequately  dealt  with  this  method  can 
be  very  accurate. 

7.1  Summary 

The  coupled  field  problem  has  been  introduced  with  the 
available  solution  methods.  Simulation  coupling  was  introduced  as 
an  branch  of  partitioned  solution  methods,  where  the  partitions  are 
created  to  decoupled  the  dynamics  into  single  fields.  This  allowed 
the  use  of  constraints  and  single  field  simulation  tools  to  solve 
coupled  field  problems. 

Integration  methods  were  examined  in  detail,  focusing  on 
linear  multistep  methods.  Operator  expressions  were  introduced  to 
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provide  a  compact  notation  for  the  analysis  tools  developed.  These 
expressions  were  applied  to  integration  methods  and  predictors  to 
find  expressions  for  the  stability  and  accuracy  of  complex  integration 
schemes. 

Numerical  methods  of  dealing  with  constraints  were  discussed, 
focusing  on  the  classical  Lagrange  mulipliers  and  the  newer 
stabilization  schemes.  These  methods  were  used  along  with  the 
integration  tools  to  develop  the  simulation  coupling  algorithm. 
Although  two  methods  of  simulation  coupling  were  introduced  the 
concurrent  evaluation  was  given  considerably  more  attention.  The 
validity  of  the  method  and  the  analysis  tools  developed  was 
examined  in  a  series  of  small  scale  problems  before  applying  the 
theory  to  the  Space  Station  simulation. 

This  theory  is  in  no  way  meant  to  be  suggested  as  superior  to 
any  other  methods  of  solving  coupled  field  problems.  It  is 
introduced  here  as  a  possible  alternative  to  preforming  large  scale 
analysis  of  complex  problems  such  as  multibody  dynamics. 

7.2  Significant  Findings  and  Results 

Simulation  coupling  has  proven  to  be  worthy  of  further 
consideration  as  an  alternative  solution  process  for  coupled  field 
problems.  It  shows  considerable  success  on  small  flexible  and  rigid 
body  problems.  It  allows  a  considerable  savings  in  time  and  effort  to 
set  up  the  simulation  for  and  find  the  solution  of  large  coupled 
problems.  It  is  not  at  the  time  a  perfect  alternative  and  needs 
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further  investigation  in  several  areas  before  application  to  larger 
problems.  However,  the  potential  shown  warrants  this  effort. 

7.3  Future  Work 

Of  the  several  areas  that  require  more  work,  the  one  with  the 
most  potential  for  aiding  the  success  of  simulation  coupling  as  a 
viable  solution  form  is  prediction  methods.  Only  the  simplest  form  of 
prediction  has  been  used  in  this  thesis  but  current  research  suggests 
that  large  gains  in  time  step  limits  and  accuracy  can  be  made 
through  more  advanced  predictors.  Work  done  in  the  past  on  other 
partitioned  solution  forms  suggests  that  ‘optimal’  predictors  may 
exist  and  that  they  have  large  impact  on  a  method’s  success. 

More  work  needs  to  be  done  in  the  area  of  handling  and 
incorporating  constraints.  This  area  is  essential  to  the  success  of 
simulation  coupling.  At  the  current  time  there  is  still  a  great  deal  of 
debate  as  to  which  of  the  many  possibilities  is  most  suitable  for 
widescale  application. 

Finally,  an  area  with  enormous  potential  for  application  of  the 
simulation  coupling  algorithm  is  parallel  computing.  Seperate 
simulation  elements  could  be  tasked  to  the  various  processors  with 
one  to  oversee  data  management  and  execution.  This  would  make 
simulation  coupling  a  potentially  cost  saving  method  not  only  in 
effort  required  to  create  simulation  codes,  but  also  in  real  time 
necessary  to  run  simulations.  Simulation  coupling  readily  lends  itself 
to  such  a  multitasking  enviroment. 
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Appendix  A 


Satellite  Rigid  Body  Plots 


As  mentioned  in  Section  5.3.3  this  Appendix  contains  all  of  the 
plots  for  the  third  simple  numerical  example.  The  particular  forcing 
conditions  for  each  of  the  four  cases  are  described  within  the 
appropriate  section.  The  third  and  fourth  cases  have  long  and  short 
duration  runs.  The  units  for  all  translations  degrees  are  meters  and 
meters  per  second,  rotations  are  degrees  and  degrees  per  second,  and 
time  is  in  seconds. 

A.l  First  Case  -  Z  Axis  Forcing 

A  constant  forcing  of  1  N  is  applied  in  the  direction  of  the  z 
axis.  There  is  an  initial  spin  rate  of  0.5  rpm  about  the  z  axis.  Due  to 
symmetry  only  a  z  axis  translation  results.  Responses  included  are 
rotor  and  platform  translational  displacements  as  well  as  the 
constraint  error  responses.  Angular  responses  are  not  included  since 
the  only  angular  displacement  that  occurs  is  a  constant  increase 
about  the  z  axis. 
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A.2  Second  Case  -  Z  Axis  Torque 

A  constant  forcing  of  10  N-m  torque  is  applied  about  the  z  axis. 
There  is  an  initial  spin  rate  of  0.5  rpm  about  the  z  axis.  Due  to 
symmetry  only  a  z  axis  rotation  results.  Responses  included  are 
rotor  and  platform  rotational  displacements  as  well  as  the  constraint 
error  responses.  Translational  responses  are  not  included  since  all 
displacements  are  zero. 


Figure  A.2-1  Rotational  Response  f 
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A.3  Third  Case  -  Tri-Axis  Forcing 

A  constant  forcing  of  1  N  is  applied  to  all  three  axes.  There  is 
an  initial  spin  rate  of  0.5  rpm  about  the  z  axis.  Due  to  symmetry 
there  is  no  coupling  into  the  rotational  DOFs  and  only  the  constant  z 
axis  rotation  results.  Responses  included  are  rotor  and  platform 
translational  displacements  and  velocities  as  well  as  the  constraint 
error  responses.  Rotational  responses  are  not  included. 

A  ten  second  run  and  a  four  minute  run  are  included.  The  ten 
second  run  displays  the  initial  damping  of  the  constraint  error.  The 
longer  run  shows  two  cycles  of  a  longer  period  error  induced  from 
the  constant  rotation  about  the  z  axis. 
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Figure  A. 3-3  Translational  Error  Response,  Short  Run,  Case  3 
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Figure  A. 3-4  Translational  Rotor  Response,  Long  Run,  Case  3 
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Velocity  Platform  ^  2  X  Velocity  Platform 


Figure  A. 3-6  Translational  Error  Response,  Long  Run,  Case  3 


123 


STRUCTURAL  SIMULATION  COUPLING  FOR  TRANSIENT  ANALYSIS 


A.4  Fourth  Case  -  Non-Spin  Axis  Torque 

A  0.1  second  torque  of  1  N-m  is  applied  about  a  non-spin  axis. 
There  is  an  initial  spin  rate  of  0.5  rpm  about  the  z  axis.  For  the  first 
time  there  is  significant  coupling  between  the  translational  and 
rotational  DOFs.  Responses  included  are  rotor  and  platform 
translational  displacements  and  velocities  as  well  as  the  rotational 
displacements  and  velocities.  Constraint  error  responses  are  also 
included. 

A  ten  second  run  and  a  sixty  second  run  are  included.  The  ten 
second  run  displays  the  initial  damping  of  the  constraint  error.  The 
longer  run  shows  a  longer  period  error  induced  from  the  rotations 
and  the  effects  of  cross-coupling.  The  error  responses  take  on  a 
slightly  different  form  than  seen  previously  due  to  the  fact  that  the 
forcing  is  cycled  so  quickly. 
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Figure  A. 4-2  Rota 
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Figure  A.4-4  Rotational  Platform  Response,  Short  Run,  Case  4 


Platform  Theta  2  Rate  Platform  Theta  1  Rate  Platform 


Error  Theta  3  g'  k  Error  Theta  2  fc  f  Error  Theta 


Z  Position  Rotor  2  £,  Y  Position  Rotor  o  a,  X  Position  Rotor 


APPENDIX 


ms 


jJTTiTF 


Ijk' 


Time  Time 

Figure  A.4-7  Translational  Rotor  Response,  Long  Run,  Case  4 
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Figure  A.4-9  Translational  Platform  Response,  Long  Run,  Case  4 
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Figure  A.4-11  Translational  Error  Response,  Long  Run,  Case  4 
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Figure  A.4-12  Rotational  Error  Response,  Long  Run,  Case  4 
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Appendix  B 


Space  Station  Plots 


As  mentioned  in  Section  7.2  this  Appendix  contains  the  plots 
for  the  Space  Station  example.  The  Appendix  is  split  into  two 
sections,  one  for  the  rigid  body  model  and  one  for  the  flexible  model. 
The  individual  cases  are  described  in  each  section.  In  all  plots  the 
dotted  line  represents  the  ‘truth’  model  response.  The  truth  model  is 
actually  the  repsonse  taken  using  a  three  body  mulitibody  solution. 

B.l  Rigid  Model  Plots 

This  section  includes  four  separate  cases,  a  non-articulated  and 
an  articulated  case  for  each  of  the  two  forcing  conditions.  The  two 
forcings  are  a  100  lb  force  in  the  direction  of  the  z  axis  and  a  2500 
ft-lb  torque  about  the  x  axis.  In  addition  to  the  four  case  responses 
there  are  plots  representing  a  typical  error  response  for  the  rigid 
cases.  The  particular  case  these  error  plots  were  taken  from  is  the 
non-articulated  z  forced  case. 
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Figure  B.l-1  Z  Force  Response  for  Rigid  Nonarticulated  Case 
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Figure  B.l-2  Z  Force  Response  for  Rigid  Nonarticulated  Case 
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Figure  B.l-3  X  Torque  Response  for  Rigid  Nonarticulated  Case 
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Figure  B.l-5  Z  Force  Response  for  Rigid  Articulated  Case 
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Z  Force  Response  for  Rigid  Articulated  Case 
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Figure  B.l-7  X  Torque  Response  for  Rigid  Articulated  Case 
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Figure  B.l-8  X  Torque  Response  for  Rigid  Articulated  Case 
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B.2  Flexible  Model  Plots 

The  flexible  model  was  run  with  10  modes  for  each  subdomain. 
The  responses  are  included  for  each  forcing  case  as  well  as  the  error 
plots  for  the  x  torque  case.  The  error  plots  for  the  z  force  case  are 
not  significantly  different  and  are  not  included. 
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Figure  B.2-3  X  Torque  Response  for  Flexible  Articulated  Case 
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Figure  B.2-4  X  Torque  Response  for  Flexible  Articulated  Case 
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Figure  B.2-5  Error  Response  for  Flexible  Nonarticulated  Case 
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Numerical  Improvements  for 
Simulation  Coupling 


Since  the  computational  dynamist  has  so  few  free  parameters 
available  to  him  some  addition  measures  are  sometimes  necessary  to 
aid  in  achieving  the  desired  stability  and  accuracy.  The  most 
common  of  these  measures  are  additional  prediction  steps  and 
iteration  loops. 

C.l  Use  of  Rigid  Body  Predictors 

One  main  difficulty  with  simulation  coupling  is  that  the 
corrective  forces  represent  a  delayed  feedback  of  information.  A 
force  applied  to  one  body  takes  one  full  time  step  before  its  effects 
can  be  felt  at  a  connecting  body.  One  easy  correction  for  coupled 
structural  problems  is  the  use  of  a  rigid  body  predictor  step. 
Consider  as  an  example  the  simple  two  field  problem  illustrated  in 
Figure  C.1-1. 
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Given  knowledge  of  the  geometry  of  the  problem  and  the 
mass/inertia  matrices  of  the  separate  fields  a  global  rigid  body  mass 
matrix,  Mg,  is  easily  assembled.  Using  this  matrix  and  the  total  force 
applied  to  the  coupled  domain,  the  rigid  body  accelerations  may  be 
found  by  solving  the  matrix  form  of  f  =  Mucm.  These  accelerations 
are  used  along  with  basic  kinematics  to  solve  for  the  accelerations 
any  point  in  the  system.  The  accelerations  at  the  boundary,  point  A, 
are 

u’a  =  Ucm  +  ojxuCm  +  coxrA  +  cox(a)xrA)  (C.l-1) 

Using  these  values  and  the  mass/inertia  matrix  of  field  A,  the  forces 
are  found  which,  when  applied  at  point  A  produce  the  accelerations 
iiA.  To  maintain  dynamic  equilibrium,  an  equal  and  opposite  set  of 
forces  must  be  applied  to  field  B.  These  forces  represent  the  set  of 
d’Alembert  forces  which  act  at  the  boundary  to  ensure  the 
constraints  are  met  under  rigid  conditions.  This  force  balance  is 
presented  in  Figure  C.l-2. 
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Figure  C.l-2  Force  Balance  Using  Rigid  Body  Prediction 


The  overall  effect  of  using  a  rigid  body  predictor  is  to 
appropriately  distribute  the  forces  applied  at  given  time  to  all  the 
coupled  fields.  In  this  manner  the  effects  of  a  force  on  a  body  are 
felt  immediately  across  the  entire  domain.  In  this  way  the  problems 
associated  with  time  delayed  information  are  avoided.  An  additional 
benefit  comes  from  the  fact  that  the  rigid  body  predictor  keeps  the 
simulation  coupled  solution  close  to  the  actual  solution.  Because  of 
this  smaller  penalties  and  larger  time  steps  than  without  the 
predictor  may  be  used. 


C.2  Example  of  Rigid  Body  Prediction 


As  an  example  using  a  rigid  body  predictor  the  Space  Station 
rigid  model  with  articulation  is  used.  The  penalties  required  to  keep 
translation  errors  less  than  10'8  meters  and  angular  errors  less  than 
10‘6  degrees  are 
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a  K 

Translation  100,000  100.0 

Rotation  50,000  10.0 

Since  the  penalties  are  so  much  lower  than  without  using  a  rigid 
body  predictor,  the  time  step  could  be  lowered  from  0.0001  to  0.005. 
The  case  is  run  using  a  100  lb  force  in  the  z  direction.  The  plots  for 
rigid  body  translations  and  rotations  are  included  as  Figures  C. 2-1 
and  C.2-2.  In  this  example  all  responses  are  basically  overlays  with 
a  small  drift  in  the  63  angle. 
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Figure  C.2-1  Z  Force  Response  for  Rigid  Articulated  Case 
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